PartBunch.cpp 44.2 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
// ------------------------------------------------------------------------
// $RCSfile: PartBunch.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class PartBunch
//   Interface to a particle bunch.
//   Can be used to avoid use of a template in user code.
//
// ------------------------------------------------------------------------
// Class category: Algorithms
// ------------------------------------------------------------------------
//
// $Date: 2004/11/12 18:57:53 $
// $Author: adelmann $
//
// ------------------------------------------------------------------------

#include "Algorithms/PartBunch.h"
#include "FixedAlgebra/FMatrix.h"
#include "FixedAlgebra/FVector.h"
#include <iostream>
#include <cfloat>
#include <fstream>
#include <iomanip>
adelmann's avatar
adelmann committed
29
#include <memory>
30
#include <utility>
gsell's avatar
gsell committed
31

32

33 34
#include "Distribution/Distribution.h"  // OPAL file
#include "Structure/FieldSolver.h"      // OPAL file
adelmann's avatar
adelmann committed
35
#include "Utilities/GeneralClassicException.h"
gsell's avatar
gsell committed
36

kraus's avatar
kraus committed
37
#include "Algorithms/ListElem.h"
gsell's avatar
gsell committed
38 39 40 41 42 43 44 45

#include <gsl/gsl_rng.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_qrng.h>

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
46 47
#include <boost/format.hpp>

48
//#define DBG_SCALARFIELD
adelmann's avatar
adelmann committed
49
//#define FIELDSTDOUT
50

gsell's avatar
gsell committed
51 52 53
// Class PartBunch
// ------------------------------------------------------------------------

54 55
PartBunch::PartBunch(const PartData *ref): // Layout is set using setSolver()
    PartBunchBase<double, 3>(new PartBunch::pbase_t(new Layout_t()), ref),
frey_m's avatar
frey_m committed
56
    interpolationCacheSet_m(false)
57
{
kraus's avatar
kraus committed
58

59
}
60

gsell's avatar
gsell committed
61

62 63 64 65
PartBunch::~PartBunch() {

}

frey_m's avatar
AMR:  
frey_m committed
66 67 68
// PartBunch::pbase_t* PartBunch::clone() {
//     return new pbase_t(new Layout_t());
// }
gsell's avatar
gsell committed
69

70

71 72 73 74
void PartBunch::initialize(FieldLayout_t *fLayout) {
    Layout_t* layout = static_cast<Layout_t*>(&getLayout());
    layout->getLayout().changeDomain(*fLayout);
}
gsell's avatar
gsell committed
75

76
void PartBunch::runTests() {
77

78 79 80
    Vector_t ll(-0.005);
    Vector_t ur(0.005);

81 82
    setBCAllPeriodic();

83
    NDIndex<3> domain = getFieldLayout().getDomain();
84
    for (unsigned int i = 0; i < Dimension; i++)
85
        nr_m[i] = domain[i].length();
86

87
    for (int i = 0; i < 3; i++)
88
        hr_m[i] = (ur[i] - ll[i]) / nr_m[i];
89

90 91
    getMesh().set_meshSpacing(&(hr_m[0]));
    getMesh().set_origin(ll);
92

93 94
    rho_m.initialize(getMesh(),
                     getFieldLayout(),
frey_m's avatar
frey_m committed
95
                     GuardCellSizes<Dimension>(1),
96 97 98
                     bc_m);
    eg_m.initialize(getMesh(),
                    getFieldLayout(),
frey_m's avatar
frey_m committed
99
                    GuardCellSizes<Dimension>(1),
100 101
                    vbc_m);

frey_m's avatar
frey_m committed
102 103 104 105 106 107
    fs_m->solver_m->test(this);
}


void PartBunch::do_binaryRepart() {
    get_bounds(rmin_m, rmax_m);
kraus's avatar
kraus committed
108

109
    pbase_t* underlyingPbase =
frey_m's avatar
frey_m committed
110
        dynamic_cast<pbase_t*>(pbase.get());
kraus's avatar
kraus committed
111

112
    BinaryRepartition(*underlyingPbase);
frey_m's avatar
frey_m committed
113 114 115
    update();
    get_bounds(rmin_m, rmax_m);
    boundp();
116 117 118
}


gsell's avatar
gsell committed
119 120 121 122 123 124 125 126 127 128 129 130
void PartBunch::computeSelfFields(int binNumber) {
    IpplTimings::startTimer(selfFieldTimer_m);

    /// Set initial charge density to zero. Create image charge
    /// potential field.
    rho_m = 0.0;
    Field_t imagePotential = rho_m;

    /// Set initial E field to zero.
    eg_m = Vector_t(0.0);

    if(fs_m->hasValidSolver()) {
131
        /// Mesh the whole domain
132
        resizeMesh();
133

gsell's avatar
gsell committed
134 135
        /// Scatter charge onto space charge grid.
        this->Q *= this->dt;
136 137
        if(!interpolationCacheSet_m) {
            if(interpolationCache_m.size() < getLocalNum()) {
138 139 140
                interpolationCache_m.create(getLocalNum() - interpolationCache_m.size());
            } else {
                interpolationCache_m.destroy(interpolationCache_m.size() - getLocalNum(),
141 142
                                             getLocalNum(),
                                             true);
143 144 145 146 147 148
            }
            interpolationCacheSet_m = true;
            this->Q.scatter(this->rho_m, this->R, IntrplCIC_t(), interpolationCache_m);
        } else {
            this->Q.scatter(this->rho_m, IntrplCIC_t(), interpolationCache_m);
        }
149

gsell's avatar
gsell committed
150 151 152 153
        this->Q /= this->dt;
        this->rho_m /= getdT();

        /// Calculate mesh-scale factor and get gamma for this specific slice (bin).
154 155
        double scaleFactor = 1;
        // double scaleFactor = Physics::c * getdT();
gsell's avatar
gsell committed
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
        double gammaz = getBinGamma(binNumber);

        /// Scale charge density to get charge density in real units. Account for
        /// Lorentz transformation in longitudinal direction.
        double tmp2 = 1 / hr_m[0] * 1 / hr_m[1] * 1 / hr_m[2] / (scaleFactor * scaleFactor * scaleFactor) / gammaz;
        rho_m *= tmp2;

        /// Scale mesh spacing to real units (meters). Lorentz transform the
        /// longitudinal direction.
        Vector_t hr_scaled = hr_m * Vector_t(scaleFactor);
        hr_scaled[2] *= gammaz;

        /// Find potential from charge in this bin (no image yet) using Poisson's
        /// equation (without coefficient: -1/(eps)).
        imagePotential = rho_m;
171

gsell's avatar
gsell committed
172
        fs_m->solver_m->computePotential(rho_m, hr_scaled);
173

gsell's avatar
gsell committed
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188
        /// Scale mesh back (to same units as particle locations.)
        rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];

        /// The scalar potential is given back in rho_m
        /// and must be converted to the right units.
        rho_m *= getCouplingConstant();

        /// IPPL Grad numerical computes gradient to find the
        /// electric field (in bin rest frame).
        eg_m = -Grad(rho_m, eg_m);

        /// Scale field. Combine Lorentz transform with conversion to proper field
        /// units.
        eg_m *= Vector_t(gammaz / (scaleFactor), gammaz / (scaleFactor), 1.0 / (scaleFactor * gammaz));

189 190 191 192 193 194 195 196 197 198 199 200
        // If desired write E-field and potential to terminal
#ifdef FIELDSTDOUT
        // Immediate debug output:
        // Output potential and e-field along the x-, y-, and z-axes
        int mx = (int)nr_m[0];
        int mx2 = (int)nr_m[0] / 2;
        int my = (int)nr_m[1];
        int my2 = (int)nr_m[1] / 2;
        int mz = (int)nr_m[2];
        int mz2 = (int)nr_m[2] / 2;

        for (int i=0; i<mx; i++ )
kraus's avatar
kraus committed
201
            *gmsg << "Bin " << binNumber
202
                  << ", Self Field along x axis E = " << eg_m[i][my2][mz2]
203 204 205
                  << ", Pot = " << rho_m[i][my2][mz2]  << endl;

        for (int i=0; i<my; i++ )
206 207
            *gmsg << "Bin " << binNumber
                  << ", Self Field along y axis E = " << eg_m[mx2][i][mz2]
208 209 210
                  << ", Pot = " << rho_m[mx2][i][mz2]  << endl;

        for (int i=0; i<mz; i++ )
211 212
            *gmsg << "Bin " << binNumber
                  << ", Self Field along z axis E = " << eg_m[mx2][my2][i]
213 214 215
                  << ", Pot = " << rho_m[mx2][my2][i]  << endl;
#endif

gsell's avatar
gsell committed
216 217 218 219
        /// Interpolate electric field at particle positions.  We reuse the
        /// cached information about where the particles are relative to the
        /// field, since the particles have not moved since this the most recent
        /// scatter operation.
220 221
        Eftmp.gather(eg_m, IntrplCIC_t(), interpolationCache_m);
        //Eftmp.gather(eg_m, this->R, IntrplCIC_t());
gsell's avatar
gsell committed
222 223 224 225 226 227 228 229

        /** Magnetic field in x and y direction induced by the electric field.
         *
         *  \f[ B_x = \gamma(B_x^{'} - \frac{beta}{c}E_y^{'}) = -\gamma \frac{beta}{c}E_y^{'} = -\frac{beta}{c}E_y \f]
         *  \f[ B_y = \gamma(B_y^{'} - \frac{beta}{c}E_x^{'}) = +\gamma \frac{beta}{c}E_x^{'} = +\frac{beta}{c}E_x \f]
         *  \f[ B_z = B_z^{'} = 0 \f]
         *
         */
230
        double betaC = std::sqrt(gammaz * gammaz - 1.0) / gammaz / Physics::c;
gsell's avatar
gsell committed
231 232 233 234 235 236 237 238 239 240

        Bf(0) = Bf(0) - betaC * Eftmp(1);
        Bf(1) = Bf(1) + betaC * Eftmp(0);

        Ef += Eftmp;

        /// Now compute field due to image charge. This is done separately as the image charge
        /// is moving to -infinity (instead of +infinity) so the Lorentz transform is different.

        /// Find z shift for shifted Green's function.
241 242 243 244
        NDIndex<3> domain = getFieldLayout().getDomain();
        Vector_t origin = rho_m.get_mesh().get_origin();
        double hz = rho_m.get_mesh().get_meshSpacing(2);
        double zshift = -(2 * origin(2) + (domain[2].first() + domain[2].last() + 1) * hz) * gammaz * scaleFactor;
gsell's avatar
gsell committed
245 246 247 248 249 250 251 252 253 254 255 256

        /// Find potential from image charge in this bin using Poisson's
        /// equation (without coefficient: -1/(eps)).
        fs_m->solver_m->computePotential(imagePotential, hr_scaled, zshift);

        /// Scale mesh back (to same units as particle locations.)
        imagePotential *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];

        /// The scalar potential is given back in rho_m
        /// and must be converted to the right units.
        imagePotential *= getCouplingConstant();

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
257
        const int dumpFreq = 100;
258
#ifdef DBG_SCALARFIELD
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
259
        VField_t tmp_eg = eg_m;
260

kraus's avatar
kraus committed
261

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
262 263 264 265 266 267 268
        if (Ippl::getNodes() == 1 && (fieldDBGStep_m + 1) % dumpFreq == 0) {
#else
        VField_t tmp_eg;

        if (false) {
#endif
            INFOMSG(level1 << "*** START DUMPING SCALAR FIELD ***" << endl);
kraus's avatar
kraus committed
269

270 271

            std::string SfileName = OpalData::getInstance()->getInputBasename();
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
272 273 274
            boost::format phi_fn("data/%1%-phi_scalar-%|2$05|.dat");
            phi_fn % SfileName % (fieldDBGStep_m / dumpFreq);

275
            std::ofstream fstr2(phi_fn.str());
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
276
            fstr2.precision(9);
277 278 279 280 281 282

            NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
            Vector_t origin = rho_m.get_mesh().get_origin();
            Vector_t spacing(rho_m.get_mesh().get_meshSpacing(0),
                             rho_m.get_mesh().get_meshSpacing(1),
                             rho_m.get_mesh().get_meshSpacing(2));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
283 284 285 286 287 288 289 290
            Vector_t rmin, rmax;
            get_bounds(rmin, rmax);

            INFOMSG(level1
                    << (rmin(0) - origin(0)) / spacing(0) << "\t"
                    << (rmin(1)  - origin(1)) / spacing(1) << "\t"
                    << (rmin(2)  - origin(2)) / spacing(2) << "\t"
                    << rmin(2) << endl);
291 292 293
            for (int x = myidx[0].first(); x <= myidx[0].last(); x++) {
                for (int y = myidx[1].first(); y <= myidx[1].last(); y++) {
                    for (int z = myidx[2].first(); z <= myidx[2].last(); z++) {
294 295 296 297 298
                        fstr2 << std::setw(5) << x + 1
                              << std::setw(5) << y + 1
                              << std::setw(5) << z + 1
                              << std::setw(17) << origin(2) + z * spacing(2)
                              << std::setw(17) << rho_m[x][y][z].get()
frey_m's avatar
frey_m committed
299
                              << std::setw(17) << imagePotential[x][y][z].get() << std::endl;
300 301 302 303 304
                    }
                }
            }
            fstr2.close();

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
305
            INFOMSG(level1 << "*** FINISHED DUMPING SCALAR FIELD ***" << endl);
306 307
        }

gsell's avatar
gsell committed
308 309 310 311 312 313 314 315 316
        /// IPPL Grad numerical computes gradient to find the
        /// electric field (in rest frame of this bin's image
        /// charge).
        eg_m = -Grad(imagePotential, eg_m);

        /// Scale field. Combine Lorentz transform with conversion to proper field
        /// units.
        eg_m *= Vector_t(gammaz / (scaleFactor), gammaz / (scaleFactor), 1.0 / (scaleFactor * gammaz));

317 318 319 320
        // If desired write E-field and potential to terminal
#ifdef FIELDSTDOUT
        // Immediate debug output:
        // Output potential and e-field along the x-, y-, and z-axes
321 322 323 324 325 326
        //int mx = (int)nr_m[0];
        //int mx2 = (int)nr_m[0] / 2;
        //int my = (int)nr_m[1];
        //int my2 = (int)nr_m[1] / 2;
        //int mz = (int)nr_m[2];
        //int mz2 = (int)nr_m[2] / 2;
327 328

        for (int i=0; i<mx; i++ )
kraus's avatar
kraus committed
329
            *gmsg << "Bin " << binNumber
330
                  << ", Image Field along x axis E = " << eg_m[i][my2][mz2]
331 332 333
                  << ", Pot = " << rho_m[i][my2][mz2]  << endl;

        for (int i=0; i<my; i++ )
334 335
            *gmsg << "Bin " << binNumber
                  << ", Image Field along y axis E = " << eg_m[mx2][i][mz2]
336 337 338
                  << ", Pot = " << rho_m[mx2][i][mz2]  << endl;

        for (int i=0; i<mz; i++ )
339 340
            *gmsg << "Bin " << binNumber
                  << ", Image Field along z axis E = " << eg_m[mx2][my2][i]
341 342 343
                  << ", Pot = " << rho_m[mx2][my2][i]  << endl;
#endif

344
#ifdef DBG_SCALARFIELD
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
345 346 347 348 349
        if (Ippl::getNodes() == 1 && (fieldDBGStep_m + 1) % dumpFreq == 0) {
#else
        if (false) {
#endif
            INFOMSG(level1 << "*** START DUMPING E FIELD ***" << endl);
350 351

            std::string SfileName = OpalData::getInstance()->getInputBasename();
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
352 353
            boost::format phi_fn("data/%1%-e_field-%|2$05|.dat");
            phi_fn % SfileName % (fieldDBGStep_m / dumpFreq);
354

355
            std::ofstream fstr2(phi_fn.str());
356 357 358 359 360 361
            fstr2.precision(9);

            Vector_t origin = eg_m.get_mesh().get_origin();
            Vector_t spacing(eg_m.get_mesh().get_meshSpacing(0),
                             eg_m.get_mesh().get_meshSpacing(1),
                             eg_m.get_mesh().get_meshSpacing(2));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
362

363
            NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
364 365 366
            for (int x = myidxx[0].first(); x <= myidxx[0].last(); x++) {
                for (int y = myidxx[1].first(); y <= myidxx[1].last(); y++) {
                    for (int z = myidxx[2].first(); z <= myidxx[2].last(); z++) {
367 368 369 370 371 372 373
                        Vector_t ef = eg_m[x][y][z].get() + tmp_eg[x][y][z].get();
                        fstr2 << std::setw(5) << x + 1
                              << std::setw(5) << y + 1
                              << std::setw(5) << z + 1
                              << std::setw(17) << origin(2) + z * spacing(2)
                              << std::setw(17) << ef(0)
                              << std::setw(17) << ef(1)
frey_m's avatar
frey_m committed
374
                              << std::setw(17) << ef(2) << std::endl;
375 376 377 378 379 380 381 382 383
                    }
                }
            }

            fstr2.close();

            //MPI_File_write_shared(file, (char*)oss.str().c_str(), oss.str().length(), MPI_CHAR, &status);
            //MPI_File_close(&file);

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
384
            INFOMSG(level1 << "*** FINISHED DUMPING E FIELD ***" << endl);
385
        }
386
        fieldDBGStep_m++;
387

gsell's avatar
gsell committed
388 389 390 391
        /// Interpolate electric field at particle positions.  We reuse the
        /// cached information about where the particles are relative to the
        /// field, since the particles have not moved since this the most recent
        /// scatter operation.
392 393
        Eftmp.gather(eg_m, IntrplCIC_t(), interpolationCache_m);
        //Eftmp.gather(eg_m, this->R, IntrplCIC_t());
gsell's avatar
gsell committed
394 395 396 397 398 399 400 401 402 403 404 405 406

        /** Magnetic field in x and y direction induced by the image charge electric field. Note that beta will have
         *  the opposite sign from the bunch charge field, as the image charge is moving in the opposite direction.
         *
         *  \f[ B_x = \gamma(B_x^{'} - \frac{beta}{c}E_y^{'}) = -\gamma \frac{beta}{c}E_y^{'} = -\frac{beta}{c}E_y \f]
         *  \f[ B_y = \gamma(B_y^{'} - \frac{beta}{c}E_x^{'}) = +\gamma \frac{beta}{c}E_x^{'} = +\frac{beta}{c}E_x \f]
         *  \f[ B_z = B_z^{'} = 0 \f]
         *
         */
        Bf(0) = Bf(0) + betaC * Eftmp(1);
        Bf(1) = Bf(1) - betaC * Eftmp(0);

        Ef += Eftmp;
407

gsell's avatar
gsell committed
408 409 410 411 412
    }
    IpplTimings::stopTimer(selfFieldTimer_m);
}

void PartBunch::resizeMesh() {
413 414 415 416
    if (fs_m->getFieldSolverType() != "SAAMG") {
        return;
    }

417 418 419 420
    double xmin = fs_m->solver_m->getXRangeMin();
    double xmax = fs_m->solver_m->getXRangeMax();
    double ymin = fs_m->solver_m->getYRangeMin();
    double ymax = fs_m->solver_m->getYRangeMax();
kraus's avatar
kraus committed
421

gsell's avatar
gsell committed
422 423 424
    if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
       ymin > rmin_m[1] || ymax < rmax_m[1]) {

425
        for (unsigned int n = 0; n < getLocalNum(); n++) {
gsell's avatar
gsell committed
426 427 428 429 430

            if(R[n](0) < xmin || R[n](0) > xmax ||
               R[n](1) < ymin || R[n](1) > ymax) {

                // delete the particle
431
                INFOMSG(level2 << "destroyed particle with id=" << ID[n] << endl;);
gsell's avatar
gsell committed
432 433 434 435 436 437 438 439 440 441
                destroy(1, n);
            }

        }

        update();
        boundp();
        get_bounds(rmin_m, rmax_m);
    }

442 443 444
    Vector_t origin = Vector_t(0.0, 0.0, 0.0);

    // update the mesh origin and mesh spacing hr_m
445
    fs_m->solver_m->resizeMesh(origin, hr_m, rmin_m, rmax_m, dh_m);
gsell's avatar
gsell committed
446 447

    getMesh().set_meshSpacing(&(hr_m[0]));
448
    getMesh().set_origin(origin);
gsell's avatar
gsell committed
449 450 451

    rho_m.initialize(getMesh(),
                     getFieldLayout(),
frey_m's avatar
frey_m committed
452
                     GuardCellSizes<Dimension>(1),
gsell's avatar
gsell committed
453 454 455
                     bc_m);
    eg_m.initialize(getMesh(),
                    getFieldLayout(),
frey_m's avatar
frey_m committed
456
                    GuardCellSizes<Dimension>(1),
gsell's avatar
gsell committed
457 458 459
                    vbc_m);

    update();
460 461

//    setGridIsFixed();
gsell's avatar
gsell committed
462 463 464 465 466 467 468 469
}

void PartBunch::computeSelfFields() {
    IpplTimings::startTimer(selfFieldTimer_m);
    rho_m = 0.0;
    eg_m = Vector_t(0.0);

    if(fs_m->hasValidSolver()) {
470
        //mesh the whole domain
471
        resizeMesh();
gsell's avatar
gsell committed
472 473 474 475 476 477 478 479 480 481

        //scatter charges onto grid
        this->Q *= this->dt;
        this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
        this->Q /= this->dt;
        this->rho_m /= getdT();

        //calculating mesh-scale factor
        double gammaz = sum(this->P)[2] / getTotalNum();
        gammaz *= gammaz;
482
        gammaz = std::sqrt(gammaz + 1.0);
483 484
        double scaleFactor = 1;
        // double scaleFactor = Physics::c * getdT();
gsell's avatar
gsell committed
485 486 487 488 489 490 491 492 493
        //and get meshspacings in real units [m]
        Vector_t hr_scaled = hr_m * Vector_t(scaleFactor);
        hr_scaled[2] *= gammaz;

        //double tmp2 = 1/hr_m[0] * 1/hr_m[1] * 1/hr_m[2] / (scaleFactor*scaleFactor*scaleFactor) / gammaz;
        double tmp2 = 1 / hr_scaled[0] * 1 / hr_scaled[1] * 1 / hr_scaled[2];
        //divide charge by a 'grid-cube' volume to get [C/m^3]
        rho_m *= tmp2;

494 495
#ifdef DBG_SCALARFIELD
        INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
frey_m's avatar
frey_m committed
496
        std::ofstream fstr1;
497 498
        fstr1.precision(9);

499
        std::string SfileName = OpalData::getInstance()->getInputBasename();
500

501
        std::string rho_fn = std::string("data/") + SfileName + std::string("-rho_scalar-") + std::to_string(fieldDBGStep_m);
frey_m's avatar
frey_m committed
502
        fstr1.open(rho_fn.c_str(), std::ios::out);
503
        NDIndex<3> myidx1 = getFieldLayout().getLocalNDIndex();
504 505 506
        for (int x = myidx1[0].first(); x <= myidx1[0].last(); x++) {
            for (int y = myidx1[1].first(); y <= myidx1[1].last(); y++) {
                for (int z = myidx1[2].first(); z <= myidx1[2].last(); z++) {
frey_m's avatar
frey_m committed
507
                    fstr1 << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  rho_m[x][y][z].get() << std::endl;
508 509 510 511 512 513
                }
            }
        }
        fstr1.close();
        INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif
514

gsell's avatar
gsell committed
515 516
        // charge density is in rho_m
        fs_m->solver_m->computePotential(rho_m, hr_scaled);
517

gsell's avatar
gsell committed
518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533
        //do the multiplication of the grid-cube volume coming
        //from the discretization of the convolution integral.
        //this is only necessary for the FFT solver
        //FIXME: later move this scaling into FFTPoissonSolver
        if(fs_m->getFieldSolverType() == "FFT" || fs_m->getFieldSolverType() == "FFTBOX")
            rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];

        // the scalar potential is given back in rho_m in units
        // [C/m] = [F*V/m] and must be divided by
        // 4*pi*\epsilon_0 [F/m] resulting in [V]
        rho_m *= getCouplingConstant();

        //write out rho
#ifdef DBG_SCALARFIELD
        INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);

frey_m's avatar
frey_m committed
534
        std::ofstream fstr2;
gsell's avatar
gsell committed
535 536
        fstr2.precision(9);

537
        std::string phi_fn = std::string("data/") + SfileName + std::string("-phi_scalar-") + std::to_string(fieldDBGStep_m);
frey_m's avatar
frey_m committed
538
        fstr2.open(phi_fn.c_str(), std::ios::out);
gsell's avatar
gsell committed
539
        NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
540 541 542
        for (int x = myidx[0].first(); x <= myidx[0].last(); x++) {
            for (int y = myidx[1].first(); y <= myidx[1].last(); y++) {
                for (int z = myidx[2].first(); z <= myidx[2].last(); z++) {
frey_m's avatar
frey_m committed
543
                    fstr2 << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  rho_m[x][y][z].get() << std::endl;
gsell's avatar
gsell committed
544 545 546 547 548 549 550 551 552 553 554 555 556 557 558
                }
            }
        }
        fstr2.close();

        INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif

        // IPPL Grad divides by hr_m [m] resulting in
        // [V/m] for the electric field
        eg_m = -Grad(rho_m, eg_m);

        eg_m *= Vector_t(gammaz / (scaleFactor), gammaz / (scaleFactor), 1.0 / (scaleFactor * gammaz));

        //write out e field
559
#ifdef FIELDSTDOUT
560 561
        // Immediate debug output:
        // Output potential and e-field along the x-, y-, and z-axes
562 563 564 565 566 567
        int mx = (int)nr_m[0];
        int mx2 = (int)nr_m[0] / 2;
        int my = (int)nr_m[1];
        int my2 = (int)nr_m[1] / 2;
        int mz = (int)nr_m[2];
        int mz2 = (int)nr_m[2] / 2;
568

569
        for (int i=0; i<mx; i++ )
570
            *gmsg << "Field along x axis Ex = " << eg_m[i][my2][mz2] << " Pot = " << rho_m[i][my2][mz2]  << endl;
571

572
        for (int i=0; i<my; i++ )
573
            *gmsg << "Field along y axis Ey = " << eg_m[mx2][i][mz2] << " Pot = " << rho_m[mx2][i][mz2]  << endl;
574

575
        for (int i=0; i<mz; i++ )
576
            *gmsg << "Field along z axis Ez = " << eg_m[mx2][my2][i] << " Pot = " << rho_m[mx2][my2][i]  << endl;
577
#endif
578

579
#ifdef DBG_SCALARFIELD
gsell's avatar
gsell committed
580 581 582 583 584
        INFOMSG("*** START DUMPING E FIELD ***" << endl);
        //ostringstream oss;
        //MPI_File file;
        //MPI_Status status;
        //MPI_Info fileinfo;
585
        //MPI_File_open(Ippl::getComm(), "rho_scalar", MPI_MODE_WRONLY | MPI_MODE_CREATE, fileinfo, &file);
frey_m's avatar
frey_m committed
586
        std::ofstream fstr;
gsell's avatar
gsell committed
587 588
        fstr.precision(9);

589
        std::string e_field = std::string("data/") + SfileName + std::string("-e_field-") + std::to_string(fieldDBGStep_m);
frey_m's avatar
frey_m committed
590
        fstr.open(e_field.c_str(), std::ios::out);
gsell's avatar
gsell committed
591
        NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
592 593 594
        for (int x = myidxx[0].first(); x <= myidxx[0].last(); x++) {
            for (int y = myidxx[1].first(); y <= myidxx[1].last(); y++) {
                for (int z = myidxx[2].first(); z <= myidxx[2].last(); z++) {
frey_m's avatar
frey_m committed
595
                    fstr << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  eg_m[x][y][z].get() << std::endl;
gsell's avatar
gsell committed
596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621
                }
            }
        }

        fstr.close();
        fieldDBGStep_m++;

        //MPI_File_write_shared(file, (char*)oss.str().c_str(), oss.str().length(), MPI_CHAR, &status);
        //MPI_File_close(&file);

        INFOMSG("*** FINISHED DUMPING E FIELD ***" << endl);
#endif

        // interpolate electric field at particle positions.  We reuse the
        // cached information about where the particles are relative to the
        // field, since the particles have not moved since this the most recent
        // scatter operation.
        Ef.gather(eg_m, this->R,  IntrplCIC_t());

        /** Magnetic field in x and y direction induced by the eletric field
         *
         *  \f[ B_x = \gamma(B_x^{'} - \frac{beta}{c}E_y^{'}) = -\gamma \frac{beta}{c}E_y^{'} = -\frac{beta}{c}E_y \f]
         *  \f[ B_y = \gamma(B_y^{'} - \frac{beta}{c}E_x^{'}) = +\gamma \frac{beta}{c}E_x^{'} = +\frac{beta}{c}E_x \f]
         *  \f[ B_z = B_z^{'} = 0 \f]
         *
         */
622
        double betaC = std::sqrt(gammaz * gammaz - 1.0) / gammaz / Physics::c;
gsell's avatar
gsell committed
623 624 625 626 627 628 629

        Bf(0) = Bf(0) - betaC * Ef(1);
        Bf(1) = Bf(1) + betaC * Ef(0);
    }
    IpplTimings::stopTimer(selfFieldTimer_m);
}

630 631 632 633 634 635 636 637 638
/**
 * \method computeSelfFields_cycl()
 * \brief Calculates the self electric field from the charge density distribution for use in cyclotrons
 * \see ParallelCyclotronTracker
 * \warning none yet
 *
 * Comments -DW:
 * I have made some changes in here:
 * -) Some refacturing to make more similar to computeSelfFields()
639
 * -) Added meanR and quaternion to be handed to the function so that SAAMG solver knows how to rotate the boundary geometry correctly.
640 641 642
 * -) Fixed an error where gamma was not taken into account correctly in direction of movement (y in cyclotron)
 * -) Comment: There is no account for image charges in the cyclotron tracker (yet?)!
 */
643
void PartBunch::computeSelfFields_cycl(double gamma) {
644

645 646 647 648 649 650 651 652
    IpplTimings::startTimer(selfFieldTimer_m);

    /// set initial charge density to zero.
    rho_m = 0.0;

    /// set initial E field to zero
    eg_m = Vector_t(0.0);

kraus's avatar
kraus committed
653
    if(fs_m->hasValidSolver()) {
654
        /// mesh the whole domain
655
        resizeMesh();
656

657 658 659 660
        /// scatter particles charge onto grid.
        this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());

        /// Lorentz transformation
661
        /// In particle rest frame, the longitudinal length (y for cyclotron) enlarged
662 663 664
        Vector_t hr_scaled = hr_m ;
        hr_scaled[1] *= gamma;

665 666 667 668 669
        /// from charge (C) to charge density (C/m^3).
        double tmp2 = 1.0 / (hr_scaled[0] * hr_scaled[1] * hr_scaled[2]);
        rho_m *= tmp2;

        // If debug flag is set, dump scalar field (charge density 'rho') into file under ./data/
670 671
#ifdef DBG_SCALARFIELD
        INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
frey_m's avatar
frey_m committed
672
        std::ofstream fstr1;
673 674
        fstr1.precision(9);

675 676 677
        std::ostringstream istr;
        istr << fieldDBGStep_m;

678
        std::string SfileName = OpalData::getInstance()->getInputBasename();
679

680
        std::string rho_fn = std::string("data/") + SfileName + std::string("-rho_scalar-") + std::string(istr.str());
frey_m's avatar
frey_m committed
681
        fstr1.open(rho_fn.c_str(), std::ios::out);
682
        NDIndex<3> myidx1 = getFieldLayout().getLocalNDIndex();
683 684 685
        for (int x = myidx1[0].first(); x <= myidx1[0].last(); x++) {
            for (int y = myidx1[1].first(); y <= myidx1[1].last(); y++) {
                for (int z = myidx1[2].first(); z <= myidx1[2].last(); z++) {
frey_m's avatar
frey_m committed
686
                    fstr1 << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  rho_m[x][y][z].get() << std::endl;
687 688 689 690 691 692
                }
            }
        }
        fstr1.close();
        INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif
693

694
        /// now charge density is in rho_m
kraus's avatar
kraus committed
695
        /// calculate Possion equation (without coefficient: -1/(eps))
696
        fs_m->solver_m->computePotential(rho_m, hr_scaled);
697 698 699 700

        //do the multiplication of the grid-cube volume coming
        //from the discretization of the convolution integral.
        //this is only necessary for the FFT solver
701
        //TODO FIXME: later move this scaling into FFTPoissonSolver
702
        if(fs_m->getFieldSolverType() == "FFT" || fs_m->getFieldSolverType() == "FFTBOX")
703
            rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];
704 705 706 707

        /// retrive coefficient: -1/(eps)
        rho_m *= getCouplingConstant();

kraus's avatar
kraus committed
708
        // If debug flag is set, dump scalar field (potential 'phi') into file under ./data/
709 710 711
#ifdef DBG_SCALARFIELD
        INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);

frey_m's avatar
frey_m committed
712
        std::ofstream fstr2;
713 714
        fstr2.precision(9);

715
        std::string phi_fn = std::string("data/") + SfileName + std::string("-phi_scalar-") + std::string(istr.str());
frey_m's avatar
frey_m committed
716
        fstr2.open(phi_fn.c_str(), std::ios::out);
717
        NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
718 719 720
        for (int x = myidx[0].first(); x <= myidx[0].last(); x++) {
            for (int y = myidx[1].first(); y <= myidx[1].last(); y++) {
                for (int z = myidx[2].first(); z <= myidx[2].last(); z++) {
frey_m's avatar
frey_m committed
721
                    fstr2 << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  rho_m[x][y][z].get() << std::endl;
722 723 724 725 726 727 728 729
                }
            }
        }
        fstr2.close();

        INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif

730 731 732
        /// calculate electric field vectors from field potential
        eg_m = -Grad(rho_m, eg_m);

733
        /// Back Lorentz transformation
kraus's avatar
kraus committed
734
        /// CAVE: y coordinate needs 1/gamma factor because IPPL function Grad() divides by
735 736 737
        /// hr_m which is not scaled appropriately with Lorentz contraction in y direction
        /// only hr_scaled is! -DW
        eg_m *= Vector_t(gamma, 1.0 / gamma, gamma);
kraus's avatar
kraus committed
738

739
#ifdef FIELDSTDOUT
740 741
        // Immediate debug output:
        // Output potential and e-field along the x-, y-, and z-axes
742 743 744 745 746 747
        int mx = (int)nr_m[0];
        int mx2 = (int)nr_m[0] / 2;
        int my = (int)nr_m[1];
        int my2 = (int)nr_m[1] / 2;
        int mz = (int)nr_m[2];
        int mz2 = (int)nr_m[2] / 2;
748

749
        for (int i=0; i<mx; i++ )
750
            *gmsg << "Field along x axis Ex = " << eg_m[i][my2][mz2] << " Pot = " << rho_m[i][my2][mz2]  << endl;
751

752
        for (int i=0; i<my; i++ )
753
            *gmsg << "Field along y axis Ey = " << eg_m[mx2][i][mz2] << " Pot = " << rho_m[mx2][i][mz2]  << endl;
754

755
        for (int i=0; i<mz; i++ )
756
            *gmsg << "Field along z axis Ez = " << eg_m[mx2][my2][i] << " Pot = " << rho_m[mx2][my2][i]  << endl;
757
#endif
758

759
#ifdef DBG_SCALARFIELD
760
        // If debug flag is set, dump vector field (electric field) into file under ./data/
761
        INFOMSG("*** START DUMPING E FIELD ***" << endl);
762 763 764 765 766
        //ostringstream oss;
        //MPI_File file;
        //MPI_Status status;
        //MPI_Info fileinfo;
        //MPI_File_open(Ippl::getComm(), "rho_scalar", MPI_MODE_WRONLY | MPI_MODE_CREATE, fileinfo, &file);
frey_m's avatar
frey_m committed
767
        std::ofstream fstr;
768
        fstr.precision(9);
769

770
        std::string e_field = std::string("data/") + SfileName + std::string("-e_field-") + std::string(istr.str());
frey_m's avatar
frey_m committed
771
        fstr.open(e_field.c_str(), std::ios::out);
772
        NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
773 774 775
        for (int x = myidxx[0].first(); x <= myidxx[0].last(); x++) {
            for (int y = myidxx[1].first(); y <= myidxx[1].last(); y++) {
                for (int z = myidxx[2].first(); z <= myidxx[2].last(); z++) {
frey_m's avatar
frey_m committed
776
                    fstr << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  eg_m[x][y][z].get() << std::endl;
777 778 779
                }
            }
        }
780 781

        fstr.close();
782 783
        fieldDBGStep_m++;

784 785 786
        //MPI_File_write_shared(file, (char*)oss.str().c_str(), oss.str().length(), MPI_CHAR, &status);
        //MPI_File_close(&file);

787 788 789
        INFOMSG("*** FINISHED DUMPING E FIELD ***" << endl);
#endif

790 791 792 793
        /// interpolate electric field at particle positions.
        Ef.gather(eg_m, this->R,  IntrplCIC_t());

        /// calculate coefficient
794 795 796
        // Relativistic E&M says gamma*v/c^2 = gamma*beta/c = sqrt(gamma*gamma-1)/c
        // but because we already transformed E_trans into the moving frame we have to
        // add 1/gamma so we are using the E_trans from the rest frame -DW
797
        double betaC = std::sqrt(gamma * gamma - 1.0) / gamma / Physics::c;
798 799 800 801 802 803

        /// calculate B field from E field
        Bf(0) =  betaC * Ef(2);
        Bf(2) = -betaC * Ef(0);

    }
kraus's avatar
kraus committed
804

805 806 807 808 809 810 811
    /*
    *gmsg << "gamma =" << gamma << endl;
    *gmsg << "dx,dy,dz =(" << hr_m[0] << ", " << hr_m[1] << ", " << hr_m[2] << ") [m] " << endl;
    *gmsg << "max of bunch is (" << rmax_m(0) << ", " << rmax_m(1) << ", " << rmax_m(2) << ") [m] " << endl;
    *gmsg << "min of bunch is (" << rmin_m(0) << ", " << rmin_m(1) << ", " << rmin_m(2) << ") [m] " << endl;
    */

812 813
    IpplTimings::stopTimer(selfFieldTimer_m);
}
gsell's avatar
gsell committed
814

815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831
/**
 * \method computeSelfFields_cycl()
 * \brief Calculates the self electric field from the charge density distribution for use in cyclotrons
 * \see ParallelCyclotronTracker
 * \warning none yet
 *
 * Overloaded version for having multiple bins with separate gamma for each bin. This is necessary
 * For multi-bunch mode.
 *
 * Comments -DW:
 * I have made some changes in here:
 * -) Some refacturing to make more similar to computeSelfFields()
 * -) Added meanR and quaternion to be handed to the function (TODO: fall back to meanR = 0 and unit quaternion
 *    if not specified) so that SAAMG solver knows how to rotate the boundary geometry correctly.
 * -) Fixed an error where gamma was not taken into account correctly in direction of movement (y in cyclotron)
 * -) Comment: There is no account for image charges in the cyclotron tracker (yet?)!
 */
832
void PartBunch::computeSelfFields_cycl(int bin) {
833 834 835 836 837 838 839 840 841 842 843 844
    IpplTimings::startTimer(selfFieldTimer_m);

    /// set initial charge dentsity to zero.
    rho_m = 0.0;

    /// set initial E field to zero
    eg_m = Vector_t(0.0);

    /// get gamma of this bin
    double gamma = getBinGamma(bin);

    if(fs_m->hasValidSolver()) {
845
        /// mesh the whole domain
846
        resizeMesh();
847 848 849 850 851

        /// scatter particles charge onto grid.
        this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());

        /// Lorentz transformation
852
        /// In particle rest frame, the longitudinal length (y for cyclotron) enlarged
853 854 855
        Vector_t hr_scaled = hr_m ;
        hr_scaled[1] *= gamma;

856 857 858 859
        /// from charge (C) to charge density (C/m^3).
        double tmp2 = 1.0 / (hr_scaled[0] * hr_scaled[1] * hr_scaled[2]);
        rho_m *= tmp2;

860 861 862
        // If debug flag is set, dump scalar field (charge density 'rho') into file under ./data/
#ifdef DBG_SCALARFIELD
        INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
frey_m's avatar
frey_m committed
863
        std::ofstream fstr1;
864 865 866 867 868
        fstr1.precision(9);

        std::ostringstream istr;
        istr << fieldDBGStep_m;

869
        std::string SfileName = OpalData::getInstance()->getInputBasename();
870

871
        std::string rho_fn = std::string("data/") + SfileName + std::string("-rho_scalar-") + std::string(istr.str());
frey_m's avatar
frey_m committed
872
        fstr1.open(rho_fn.c_str(), std::ios::out);
873
        NDIndex<3> myidx1 = getFieldLayout().getLocalNDIndex();
874 875 876
        for (int x = myidx1[0].first(); x <= myidx1[0].last(); x++) {
            for (int y = myidx1[1].first(); y <= myidx1[1].last(); y++) {
                for (int z = myidx1[2].first(); z <= myidx1[2].last(); z++) {
frey_m's avatar
frey_m committed
877
                    fstr1 << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  rho_m[x][y][z].get() << std::endl;
878 879 880 881 882 883 884
                }
            }
        }
        fstr1.close();
        INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif

885 886 887 888
        /// now charge density is in rho_m
        /// calculate Possion equation (without coefficient: -1/(eps))
        fs_m->solver_m->computePotential(rho_m, hr_scaled);

889 890
        // Do the multiplication of the grid-cube volume coming from the discretization of the convolution integral.
        // This is only necessary for the FFT solver. FIXME: later move this scaling into FFTPoissonSolver
891 892
        if(fs_m->getFieldSolverType() == "FFT" || fs_m->getFieldSolverType() == "FFTBOX")
            rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];
893 894 895 896

        /// retrive coefficient: -1/(eps)
        rho_m *= getCouplingConstant();

kraus's avatar
kraus committed
897
        // If debug flag is set, dump scalar field (potential 'phi') into file under ./data/
898 899 900
#ifdef DBG_SCALARFIELD
        INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);

frey_m's avatar
frey_m committed
901
        std::ofstream fstr2;
902 903
        fstr2.precision(9);

904
        std::string phi_fn = std::string("data/") + SfileName + std::string("-phi_scalar-") + std::string(istr.str());
frey_m's avatar
frey_m committed
905
        fstr2.open(phi_fn.c_str(), std::ios::out);
906
        NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
907 908 909
        for (int x = myidx[0].first(); x <= myidx[0].last(); x++) {
            for (int y = myidx[1].first(); y <= myidx[1].last(); y++) {
                for (int z = myidx[2].first(); z <= myidx[2].last(); z++) {
frey_m's avatar
frey_m committed
910
                    fstr2 << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  rho_m[x][y][z].get() << std::endl;
911 912 913 914 915 916 917 918
                }
            }
        }
        fstr2.close();

        INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif

919 920 921
        /// calculate electric field vectors from field potential
        eg_m = -Grad(rho_m, eg_m);

922
        /// Back Lorentz transformation
kraus's avatar
kraus committed
923
        /// CAVE: y coordinate needs 1/gamma factor because IPPL function Grad() divides by
924 925 926
        /// hr_m which is not scaled appropriately with Lorentz contraction in y direction
        /// only hr_scaled is! -DW
        eg_m *= Vector_t(gamma, 1.0 / gamma, gamma);
927

928
#ifdef FIELDSTDOUT
929 930
        // Immediate debug output:
        // Output potential and e-field along the x-, y-, and z-axes
931 932 933 934 935 936 937 938
        int mx = (int)nr_m[0];
        int mx2 = (int)nr_m[0] / 2;
        int my = (int)nr_m[1];
        int my2 = (int)nr_m[1] / 2;
        int mz = (int)nr_m[2];
        int mz2 = (int)nr_m[2] / 2;

        for (int i=0; i<mx; i++ )
kraus's avatar
kraus committed
939
            *gmsg << "Bin " << bin
940
                  << ", Field along x axis Ex = " << eg_m[i][my2][mz2]
941 942 943
                  << ", Pot = " << rho_m[i][my2][mz2]  << endl;

        for (int i=0; i<my; i++ )
944 945
            *gmsg << "Bin " << bin
                  << ", Field along y axis Ey = " << eg_m[mx2][i][mz2]
946 947 948
                  << ", Pot = " << rho_m[mx2][i][mz2]  << endl;

        for (int i=0; i<mz; i++ )
949 950
            *gmsg << "Bin " << bin
                  << ", Field along z axis Ez = " << eg_m[mx2][my2][i]
951 952
                  << ", Pot = " << rho_m[mx2][my2][i]  << endl;
#endif
953