PartBunch.cpp 83.8 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 "AbstractObjects/OpalData.h"   // OPAL file
#include "Distribution/Distribution.h"  // OPAL file
#include "Structure/FieldSolver.h"      // OPAL file
kraus's avatar
kraus committed
35 36
#include "Structure/LossDataSink.h"
#include "Utilities/Options.h"
adelmann's avatar
adelmann committed
37
#include "Utilities/GeneralClassicException.h"
38
#include "Utilities/Util.h"
gsell's avatar
gsell committed
39

kraus's avatar
kraus committed
40
#include "Algorithms/ListElem.h"
gsell's avatar
gsell committed
41 42 43 44 45 46 47 48

#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
49 50
#include <boost/format.hpp>

51
//#define DBG_SCALARFIELD
adelmann's avatar
adelmann committed
52
//#define FIELDSTDOUT
53

gsell's avatar
gsell committed
54 55
using Physics::pi;

56 57
using namespace std;

gsell's avatar
gsell committed
58 59 60 61 62 63 64 65 66
extern Inform *gmsg;

// Class PartBunch
// ------------------------------------------------------------------------

PartBunch::PartBunch(const PartData *ref):
    myNode_m(Ippl::myNode()),
    nodes_m(Ippl::getNodes()),
    fixed_grid(false),
67
    pbin_m(nullptr),
68 69 70
    lossDs_m(nullptr),
    pmsg_m(nullptr),
    f_stream(nullptr),
71
    lowParticleCount_m(false),
gsell's avatar
gsell committed
72
    reference(ref),
73 74
    unit_state_(units),
    stateOfLastBoundP_(unitless),
gsell's avatar
gsell committed
75 76 77 78 79
    moments_m(),
    dt_m(0.0),
    t_m(0.0),
    eKin_m(0.0),
    dE_m(0.0),
80
    spos_m(0.0),
81 82
    globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
    globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
gsell's avatar
gsell committed
83 84 85 86 87 88 89 90 91 92 93 94 95
    rmax_m(0.0),
    rmin_m(0.0),
    rrms_m(0.0),
    prms_m(0.0),
    rmean_m(0.0),
    pmean_m(0.0),
    eps_m(0.0),
    eps_norm_m(0.0),
    rprms_m(0.0),
    Dx_m(0.0),
    Dy_m(0.0),
    DDx_m(0.0),
    DDy_m(0.0),
96
    hr_m(-1.0),
gsell's avatar
gsell committed
97
    nr_m(0),
98
    fs_m(nullptr),
gsell's avatar
gsell committed
99 100
    couplingConstant_m(0.0),
    qi_m(0.0),
101
    interpolationCacheSet_m(false),
gsell's avatar
gsell committed
102 103
    distDump_m(0),
    fieldDBGStep_m(0),
104
    dh_m(1e-12),
gsell's avatar
gsell committed
105
    tEmission_m(0.0),
106 107
    bingamma_m(nullptr),
    binemitted_m(nullptr),
gsell's avatar
gsell committed
108 109
    lPath_m(0.0),
    stepsPerTurn_m(0),
110 111
    localTrackStep_m(0),
    globalTrackStep_m(0),
gsell's avatar
gsell committed
112 113
    numBunch_m(1),
    SteptoLastInj_m(0),
114 115
    partPerNode_m(nullptr),
    globalPartPerNode_m(nullptr),
116
    minLocNum_m(0),
117
    dist_m(nullptr),
118
    dcBeam_m(false) {
gsell's avatar
gsell committed
119 120 121
    addAttribute(P);
    addAttribute(Q);
    addAttribute(M);
122
    addAttribute(Phi);
gsell's avatar
gsell committed
123 124 125 126 127 128 129 130 131 132
    addAttribute(Ef);
    addAttribute(Eftmp);

    addAttribute(Bf);
    addAttribute(Bin);
    addAttribute(dt);
    addAttribute(PType);
    addAttribute(TriID);

    boundpTimer_m = IpplTimings::getTimer("Boundingbox");
133
    statParamTimer_m = IpplTimings::getTimer("Compute Statistics");
kraus's avatar
kraus committed
134 135
    selfFieldTimer_m = IpplTimings::getTimer("SelfField total");
    compPotenTimer_m  = IpplTimings::getTimer("SF: Potential");
gsell's avatar
gsell committed
136 137 138

    histoTimer_m = IpplTimings::getTimer("Histogram");

kraus's avatar
kraus committed
139 140
    distrCreate_m = IpplTimings::getTimer("Create Distr");
    distrReload_m = IpplTimings::getTimer("Load Distr");
gsell's avatar
gsell committed
141 142


143 144
    partPerNode_m = std::unique_ptr<size_t[]>(new size_t[Ippl::getNodes()]);
    globalPartPerNode_m = std::unique_ptr<size_t[]>(new size_t[Ippl::getNodes()]);
gsell's avatar
gsell committed
145

adelmann's avatar
adelmann committed
146
    lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(std::string("GlobalLosses"), !Options::asciidump));
gsell's avatar
gsell committed
147

148
    pmsg_m.release();
149 150
    //    f_stream.release();
    /*
151 152 153 154 155
      if(Ippl::getNodes() == 1) {
          f_stream = std::unique_ptr<ofstream>(new ofstream);
          f_stream->open("data/dist.dat", ios::out);
          pmsg_m = std::unique_ptr<Inform>(new Inform(0, *f_stream, 0));
      }
156
    */
157 158

    // set the default IPPL behaviour
159
    setMinimumNumberOfParticlesPerCore(0);
gsell's avatar
gsell committed
160 161 162
}

PartBunch::PartBunch(const PartBunch &rhs):
163
    IpplParticleBase<ParticleSpatialLayout<double, 3u> >(rhs),
gsell's avatar
gsell committed
164 165 166
    myNode_m(Ippl::myNode()),
    nodes_m(Ippl::getNodes()),
    fixed_grid(rhs.fixed_grid),
167
    pbin_m(nullptr),
168 169 170
    lossDs_m(nullptr),
    pmsg_m(nullptr),
    f_stream(nullptr),
171
    lowParticleCount_m(rhs.lowParticleCount_m),
gsell's avatar
gsell committed
172
    reference(rhs.reference),
173 174
    unit_state_(rhs.unit_state_),
    stateOfLastBoundP_(rhs.stateOfLastBoundP_),
gsell's avatar
gsell committed
175 176 177 178 179
    moments_m(rhs.moments_m),
    dt_m(rhs.dt_m),
    t_m(rhs.t_m),
    eKin_m(rhs.eKin_m),
    dE_m(rhs.dE_m),
180
    spos_m(0.0),
181 182
    globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
    globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
gsell's avatar
gsell committed
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
    rmax_m(rhs.rmax_m),
    rmin_m(rhs.rmin_m),
    rrms_m(rhs.rrms_m),
    prms_m(rhs.prms_m),
    rmean_m(rhs.rmean_m),
    pmean_m(rhs.pmean_m),
    eps_m(rhs.eps_m),
    eps_norm_m(rhs.eps_norm_m),
    rprms_m(rhs.rprms_m),
    Dx_m(rhs.Dx_m),
    Dy_m(rhs.Dy_m),
    DDx_m(rhs.DDx_m),
    DDy_m(rhs.DDy_m),
    hr_m(rhs.hr_m),
    nr_m(rhs.nr_m),
198
    fs_m(nullptr),
gsell's avatar
gsell committed
199 200
    couplingConstant_m(rhs.couplingConstant_m),
    qi_m(rhs.qi_m),
201
    interpolationCacheSet_m(rhs.interpolationCacheSet_m),
gsell's avatar
gsell committed
202 203 204 205
    distDump_m(rhs.distDump_m),
    fieldDBGStep_m(rhs.fieldDBGStep_m),
    dh_m(rhs.dh_m),
    tEmission_m(rhs.tEmission_m),
206 207
    bingamma_m(nullptr),
    binemitted_m(nullptr),
gsell's avatar
gsell committed
208 209
    lPath_m(rhs.lPath_m),
    stepsPerTurn_m(rhs.stepsPerTurn_m),
210 211
    localTrackStep_m(rhs.localTrackStep_m),
    globalTrackStep_m(rhs.globalTrackStep_m),
gsell's avatar
gsell committed
212 213
    numBunch_m(rhs.numBunch_m),
    SteptoLastInj_m(rhs.SteptoLastInj_m),
214 215
    partPerNode_m(nullptr),
    globalPartPerNode_m(nullptr),
216
    minLocNum_m(rhs.minLocNum_m),
217
    dist_m(nullptr),
218
    dcBeam_m(rhs.dcBeam_m) {
gsell's avatar
gsell committed
219
    ERRORMSG("should not be here: PartBunch::PartBunch(const PartBunch &rhs):" << endl);
220
    std::exit(0);
gsell's avatar
gsell committed
221 222 223 224 225 226
}

PartBunch::PartBunch(const std::vector<Particle> &rhs, const PartData *ref):
    myNode_m(Ippl::myNode()),
    nodes_m(Ippl::getNodes()),
    fixed_grid(false),
227
    pbin_m(nullptr),
228 229 230
    lossDs_m(nullptr),
    pmsg_m(nullptr),
    f_stream(nullptr),
231
    lowParticleCount_m(false),
gsell's avatar
gsell committed
232
    reference(ref),
233 234
    unit_state_(units),
    stateOfLastBoundP_(unitless),
gsell's avatar
gsell committed
235 236 237 238 239
    moments_m(),
    dt_m(0.0),
    t_m(0.0),
    eKin_m(0.0),
    dE_m(0.0),
240
    spos_m(0.0),
241 242
    globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
    globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
gsell's avatar
gsell committed
243 244 245 246 247 248 249 250 251 252 253 254 255
    rmax_m(0.0),
    rmin_m(0.0),
    rrms_m(0.0),
    prms_m(0.0),
    rmean_m(0.0),
    pmean_m(0.0),
    eps_m(0.0),
    eps_norm_m(0.0),
    rprms_m(0.0),
    Dx_m(0.0),
    Dy_m(0.0),
    DDx_m(0.0),
    DDy_m(0.0),
256
    hr_m(-1.0),
gsell's avatar
gsell committed
257
    nr_m(0),
258
    fs_m(nullptr),
gsell's avatar
gsell committed
259 260
    couplingConstant_m(0.0),
    qi_m(0.0),
261
    interpolationCacheSet_m(false),
gsell's avatar
gsell committed
262 263
    distDump_m(0),
    fieldDBGStep_m(0),
264
    dh_m(1e-12),
gsell's avatar
gsell committed
265
    tEmission_m(0.0),
266 267
    bingamma_m(nullptr),
    binemitted_m(nullptr),
gsell's avatar
gsell committed
268 269
    lPath_m(0.0),
    stepsPerTurn_m(0),
270 271
    localTrackStep_m(0),
    globalTrackStep_m(0),
gsell's avatar
gsell committed
272 273
    numBunch_m(1),
    SteptoLastInj_m(0),
274 275
    partPerNode_m(nullptr),
    globalPartPerNode_m(nullptr),
276
    minLocNum_m(0),
277
    dist_m(nullptr),
278
    dcBeam_m(false) {
gsell's avatar
gsell committed
279 280 281
    ERRORMSG("should not be here: PartBunch::PartBunch(const std::vector<Particle> &rhs, const PartData *ref):" << endl);
}

282 283 284 285
PartBunch::~PartBunch() {

}

gsell's avatar
gsell committed
286 287 288
/// \brief Need Ek for the Schottky effect calculation (eV)
double PartBunch::getEkin() const {
    if(dist_m)
289
        return dist_m->getEkin();
gsell's avatar
gsell committed
290 291 292 293 294 295 296
    else
        return 0.0;
}

/// \brief Need the work function for the Schottky effect calculation (eV)
double PartBunch::getWorkFunctionRf() const {
    if(dist_m)
297
        return dist_m->getWorkFunctionRf();
gsell's avatar
gsell committed
298 299 300 301 302 303
    else
        return 0.0;
}
/// \brief Need the laser energy for the Schottky effect calculation (eV)
double PartBunch::getLaserEnergy() const {
    if(dist_m)
304
        return dist_m->getLaserEnergy();
gsell's avatar
gsell committed
305 306 307 308
    else
        return 0.0;
}

309 310 311
/// \brief Return the fieldsolver type if we have a fieldsolver
std::string PartBunch::getFieldSolverType() const {
    if(fs_m)
kraus's avatar
kraus committed
312
        return fs_m->getFieldSolverType();
313 314 315
    else
        return "";
}
gsell's avatar
gsell committed
316

317
void PartBunch::runTests() {
318

319 320 321
    Vector_t ll(-0.005);
    Vector_t ur(0.005);

322 323
    setBCAllPeriodic();

324 325 326
    NDIndex<3> domain = getFieldLayout().getDomain();
    for(int i = 0; i < Dim; i++)
        nr_m[i] = domain[i].length();
327

328 329
    for(int i = 0; i < 3; i++)
        hr_m[i] = (ur[i] - ll[i]) / nr_m[i];
330

331 332
    getMesh().set_meshSpacing(&(hr_m[0]));
    getMesh().set_origin(ll);
333

334 335 336 337 338 339 340 341 342 343 344 345 346
    rho_m.initialize(getMesh(),
                     getFieldLayout(),
                     GuardCellSizes<Dim>(1),
                     bc_m);
    eg_m.initialize(getMesh(),
                    getFieldLayout(),
                    GuardCellSizes<Dim>(1),
                    vbc_m);

    fs_m->solver_m->test(*this);
}


gsell's avatar
gsell committed
347 348 349 350 351 352 353 354 355 356
/** \brief After each Schottky scan we delete all the particles.

 */
void PartBunch::cleanUpParticles() {

    size_t np = getTotalNum();
    bool scan = false;

    destroy(getLocalNum(), 0, true);

357
    dist_m->createOpalT(*this, np, scan);
gsell's avatar
gsell committed
358 359 360 361

    update();
}

362 363 364 365 366
void PartBunch::setDistribution(Distribution *d,
                                std::vector<Distribution *> addedDistributions,
                                size_t &np,
                                bool scan) {
    Inform m("setDistribution " );
gsell's avatar
gsell committed
367
    dist_m = d;
368

369
    dist_m->createOpalT(*this, addedDistributions, np, scan);
370

371
//    if (Options::cZero)
372
//        dist_m->create(*this, addedDistributions, np / 2, scan);
373
//    else
374
//        dist_m->create(*this, addedDistributions, np, scan);
gsell's avatar
gsell committed
375 376 377 378 379 380 381 382 383
}

void PartBunch::resetIfScan()
/*
  In case of a scan we have
  to reset some variables
 */
{
    dt = 0.0;
384
    localTrackStep_m = 0;
gsell's avatar
gsell committed
385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432
}



bool PartBunch::hasFieldSolver() {
    if(fs_m)
        return fs_m->hasValidSolver();
    else
        return false;
}

double PartBunch::getPx(int i) {
    return 0.0;
}

double PartBunch::getPy(int i) {
    return 0.0;
}

double PartBunch::getPz(int i) {
    return 0.0;
}

//ff
double PartBunch::getX(int i) {
    return this->R[i](0);
}

//ff
double PartBunch::getY(int i) {
    return this->R[i](1);
}

//ff
double PartBunch::getX0(int i) {
    return 0.0;
}

//ff
double PartBunch::getY0(int i) {
    return 0.0;
}

//ff
double PartBunch::getZ(int i) {
    return this->R[i](2);
}

433 434 435 436 437 438 439 440 441 442

/**
 * \method calcLineDensity()
 * \brief calculates the 1d line density (not normalized) and append it to a file.
 * \see ParallelTTracker
 * \warning none yet
 *
 * DETAILED TODO
 *
 */
443 444 445 446 447 448 449
void PartBunch::calcLineDensity(unsigned int nBins, std::vector<double> &lineDensity, std::pair<double, double> &meshInfo) {
    Vector_t rmin, rmax;
    get_bounds(rmin, rmax);

    if (nBins < 2) {
        NDIndex<3> grid = getFieldLayout().getDomain();
        nBins = grid[2].length();
gsell's avatar
gsell committed
450 451
    }

452 453 454 455 456
    double length = rmax(2) - rmin(2);
    double zmin = rmin(2) - dh_m * length, zmax = rmax(2) + dh_m * length;
    double hz = (zmax - zmin) / (nBins - 2);
    double perMeter = 1.0 / hz;//(zmax - zmin);
    zmin -= hz;
gsell's avatar
gsell committed
457

458 459
    lineDensity.resize(nBins, 0.0);
    std::fill(lineDensity.begin(), lineDensity.end(), 0.0);
gsell's avatar
gsell committed
460

461 462 463 464 465
    const unsigned int lN = getLocalNum();
    for (unsigned int i = 0; i < lN; ++ i) {
        const double z = R[i](2) - 0.5 * hz;
        unsigned int idx = (z - zmin) / hz;
        double tau = z - (zmin + idx * hz);
gsell's avatar
gsell committed
466

467 468
        lineDensity[idx] += Q[i] * (1.0 - tau) * perMeter;
        lineDensity[idx + 1] += Q[i] * tau * perMeter;
gsell's avatar
gsell committed
469 470
    }

471 472 473 474
    reduce(&(lineDensity[0]), &(lineDensity[0]) + nBins, &(lineDensity[0]), OpAddAssign());

    meshInfo.first = zmin;
    meshInfo.second = hz;
gsell's avatar
gsell committed
475 476 477 478
}

void PartBunch::calcGammas() {

479
    const int emittedBins = dist_m->getNumberOfEnergyBins();
480

gsell's avatar
gsell committed
481 482 483 484 485 486 487 488
    size_t s = 0;

    for(int i = 0; i < emittedBins; i++)
        bingamma_m[i] = 0.0;

    for(unsigned int n = 0; n < getLocalNum(); n++)
        bingamma_m[this->Bin[n]] += sqrt(1.0 + dot(this->P[n], this->P[n]));

489 490 491
    std::unique_ptr<size_t[]> particlesInBin(new size_t[emittedBins]);
    reduce(bingamma_m.get(), bingamma_m.get() + emittedBins, bingamma_m.get(), OpAddAssign());
    reduce(binemitted_m.get(), binemitted_m.get() + emittedBins, particlesInBin.get(), OpAddAssign());
gsell's avatar
gsell committed
492
    for(int i = 0; i < emittedBins; i++) {
493
        size_t &pInBin = particlesInBin[i];
gsell's avatar
gsell committed
494 495
        if(pInBin != 0) {
            bingamma_m[i] /= pInBin;
496
            INFOMSG(level2 << "Bin " << std::setw(3) << i << " gamma = " << std::setw(8) << std::scientific << std::setprecision(5) << bingamma_m[i] << "; NpInBin= " << std::setw(8) << std::setfill(' ') << pInBin << endl);
gsell's avatar
gsell committed
497 498
        } else {
            bingamma_m[i] = 1.0;
499
            INFOMSG(level2 << "Bin " << std::setw(3) << i << " has no particles " << endl);
gsell's avatar
gsell committed
500 501 502
        }
        s += pInBin;
    }
503
    particlesInBin.reset();
kraus's avatar
kraus committed
504 505


506
    if(s != getTotalNum() && !OpalData::getInstance()->hasGlobalGeometry())
gsell's avatar
gsell committed
507 508 509 510 511
        ERRORMSG("sum(Bins)= " << s << " != sum(R)= " << getTotalNum() << endl;);

    if(emittedBins >= 2) {
        for(int i = 1; i < emittedBins; i++) {
            if(binemitted_m[i - 1] != 0 && binemitted_m[i] != 0)
kraus's avatar
kraus committed
512 513
                INFOMSG(level2 << "d(gamma)= " << 100.0 * std::abs(bingamma_m[i - 1] - bingamma_m[i]) / bingamma_m[i] << " [%] "
                        << "between bin " << i - 1 << " and " << i << endl);
gsell's avatar
gsell committed
514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532
        }
    }
}


void PartBunch::calcGammas_cycl() {

    const int emittedBins = pbin_m->getLastemittedBin();

    for(int i = 0; i < emittedBins; i++)
        bingamma_m[i] = 0.0;
    for(unsigned int n = 0; n < getLocalNum(); n++)
        bingamma_m[this->Bin[n]] += sqrt(1.0 + dot(this->P[n], this->P[n]));
    for(int i = 0; i < emittedBins; i++) {
        reduce(bingamma_m[i], bingamma_m[i], OpAddAssign());
        if(pbin_m->getTotalNumPerBin(i) > 0)
            bingamma_m[i] /= pbin_m->getTotalNumPerBin(i);
        else
            bingamma_m[i] = 0.0;
533
        INFOMSG("Bin " << i << " : particle number = " << pbin_m->getTotalNumPerBin(i) << " gamma = " << bingamma_m[i] << endl);
gsell's avatar
gsell committed
534 535 536 537
    }

}

538
size_t PartBunch::calcNumPartsOutside(Vector_t x) {
kraus's avatar
kraus committed
539

540
    partPerNode_m[Ippl::myNode()] = 0;
kraus's avatar
kraus committed
541 542
    const Vector_t meanR = get_rmean();

543
    for(unsigned long k = 0; k < getLocalNum(); ++ k)
kraus's avatar
kraus committed
544 545 546 547 548 549 550
        if (abs(R[k](0) - meanR(0)) > x(0) ||
            abs(R[k](1) - meanR(1)) > x(1) ||
            abs(R[k](2) - meanR(2)) > x(2)) {

            ++ partPerNode_m[Ippl::myNode()];
        }

551 552
    reduce(partPerNode_m.get(), partPerNode_m.get() + Ippl::getNodes(), globalPartPerNode_m.get(), OpAddAssign());

kraus's avatar
kraus committed
553
    return *globalPartPerNode_m.get();
554
}
gsell's avatar
gsell committed
555 556 557 558 559 560 561 562 563 564 565 566 567

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()) {
568 569 570
        /// Mesh the whole domain
        if(fs_m->getFieldSolverType() == "SAAMG")
            resizeMesh();
571

gsell's avatar
gsell committed
572 573
        /// Scatter charge onto space charge grid.
        this->Q *= this->dt;
574 575
        if(!interpolationCacheSet_m) {
            if(interpolationCache_m.size() < getLocalNum()) {
576 577 578
                interpolationCache_m.create(getLocalNum() - interpolationCache_m.size());
            } else {
                interpolationCache_m.destroy(interpolationCache_m.size() - getLocalNum(),
579 580
                                             getLocalNum(),
                                             true);
581 582 583 584 585 586
            }
            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);
        }
587

gsell's avatar
gsell committed
588 589 590 591
        this->Q /= this->dt;
        this->rho_m /= getdT();

        /// Calculate mesh-scale factor and get gamma for this specific slice (bin).
592 593
        double scaleFactor = 1;
        // double scaleFactor = Physics::c * getdT();
gsell's avatar
gsell committed
594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609
        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)).
        IpplTimings::startTimer(compPotenTimer_m);
        imagePotential = rho_m;
610

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

gsell's avatar
gsell committed
613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629
        IpplTimings::stopTimer(compPotenTimer_m);

        /// 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));

630 631 632 633 634 635 636 637 638 639 640 641
        // 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++ )
642 643
	    *gmsg << "Bin " << binNumber
                  << ", Self Field along x axis E = " << eg_m[i][my2][mz2]
644 645 646
                  << ", Pot = " << rho_m[i][my2][mz2]  << endl;

        for (int i=0; i<my; i++ )
647 648
            *gmsg << "Bin " << binNumber
                  << ", Self Field along y axis E = " << eg_m[mx2][i][mz2]
649 650 651
                  << ", Pot = " << rho_m[mx2][i][mz2]  << endl;

        for (int i=0; i<mz; i++ )
652 653
            *gmsg << "Bin " << binNumber
                  << ", Self Field along z axis E = " << eg_m[mx2][my2][i]
654 655 656
                  << ", Pot = " << rho_m[mx2][my2][i]  << endl;
#endif

gsell's avatar
gsell committed
657 658 659 660
        /// 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.
661 662
        Eftmp.gather(eg_m, IntrplCIC_t(), interpolationCache_m);
        //Eftmp.gather(eg_m, this->R, IntrplCIC_t());
gsell's avatar
gsell committed
663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681

        /** 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]
         *
         */
        double betaC = sqrt(gammaz * gammaz - 1.0) / gammaz / Physics::c;

        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.
682 683 684 685
        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
686 687 688 689 690 691 692 693 694 695 696 697 698 699

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

        /// 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
700
        const int dumpFreq = 100;
701
#ifdef DBG_SCALARFIELD
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
702
        VField_t tmp_eg = eg_m;
703

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
704 705 706 707 708 709 710
        if (Ippl::getNodes() == 1 && (fieldDBGStep_m + 1) % dumpFreq == 0) {
#else
        VField_t tmp_eg;

        if (false) {
#endif
            INFOMSG(level1 << "*** START DUMPING SCALAR FIELD ***" << endl);
711 712

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

            ofstream fstr2(phi_fn.str());
            fstr2.precision(9);
718 719 720 721 722 723

            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
724 725 726 727 728 729 730 731
            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);
732 733 734 735 736 737 738 739 740 741 742 743 744 745
            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++) {
                        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()
                              << std::setw(17) << imagePotential[x][y][z].get() << endl;
                    }
                }
            }
            fstr2.close();

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
746
            INFOMSG(level1 << "*** FINISHED DUMPING SCALAR FIELD ***" << endl);
747 748
        }

gsell's avatar
gsell committed
749 750 751 752 753 754 755 756 757
        /// 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));

758 759 760 761
        // 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
762 763 764 765 766 767
        //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;
768 769

        for (int i=0; i<mx; i++ )
770 771
	    *gmsg << "Bin " << binNumber
                  << ", Image Field along x axis E = " << eg_m[i][my2][mz2]
772 773 774
                  << ", Pot = " << rho_m[i][my2][mz2]  << endl;

        for (int i=0; i<my; i++ )
775 776
            *gmsg << "Bin " << binNumber
                  << ", Image Field along y axis E = " << eg_m[mx2][i][mz2]
777 778 779
                  << ", Pot = " << rho_m[mx2][i][mz2]  << endl;

        for (int i=0; i<mz; i++ )
780 781
            *gmsg << "Bin " << binNumber
                  << ", Image Field along z axis E = " << eg_m[mx2][my2][i]
782 783 784
                  << ", Pot = " << rho_m[mx2][my2][i]  << endl;
#endif

785
#ifdef DBG_SCALARFIELD
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
786 787 788 789 790
        if (Ippl::getNodes() == 1 && (fieldDBGStep_m + 1) % dumpFreq == 0) {
#else
        if (false) {
#endif
            INFOMSG(level1 << "*** START DUMPING E FIELD ***" << endl);
791 792

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
796
            ofstream fstr2(phi_fn.str());
797 798 799 800 801 802
            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
803

804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824
            NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
            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++) {
                        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)
                              << std::setw(17) << ef(2) << endl;
                    }
                }
            }

            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
825
            INFOMSG(level1 << "*** FINISHED DUMPING E FIELD ***" << endl);
826
        }
827
        fieldDBGStep_m++;
828

gsell's avatar
gsell committed
829 830 831 832
        /// 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.
833 834
        Eftmp.gather(eg_m, IntrplCIC_t(), interpolationCache_m);
        //Eftmp.gather(eg_m, this->R, IntrplCIC_t());
gsell's avatar
gsell committed
835 836 837 838 839 840 841 842 843 844 845 846 847

        /** 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;
848

gsell's avatar
gsell committed
849 850 851 852 853
    }
    IpplTimings::stopTimer(selfFieldTimer_m);
}

void PartBunch::resizeMesh() {
854 855 856 857
    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();
858 859
    double zmin = fs_m->solver_m->getZRangeMin();
    double zmax = fs_m->solver_m->getZRangeMax();
kraus's avatar
kraus committed
860

gsell's avatar
gsell committed
861 862 863 864 865 866 867 868 869
    if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
       ymin > rmin_m[1] || ymax < rmax_m[1]) {

        for(unsigned int n = 0; n < getLocalNum(); n++) {

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

                // delete the particle
870
                INFOMSG(level2 << "destroyed particle with id=" << ID[n] << endl;);
gsell's avatar
gsell committed
871 872 873 874 875 876 877 878 879
                destroy(1, n);
            }

        }

        update();
        boundp();
        get_bounds(rmin_m, rmax_m);
    }
880 881
    Vector_t mymin = Vector_t(xmin, ymin , zmin);
    Vector_t mymax = Vector_t(xmax, ymax , zmax);
gsell's avatar
gsell committed
882

kraus's avatar
kraus committed
883
    for(int i = 0; i < 3; i++)
884
        hr_m[i]   = (mymax[i] - mymin[i])/nr_m[i];
gsell's avatar
gsell committed
885 886

    getMesh().set_meshSpacing(&(hr_m[0]));
kraus's avatar
kraus committed
887
    getMesh().set_origin(mymin);
gsell's avatar
gsell committed
888 889 890 891 892 893 894 895 896 897 898

    rho_m.initialize(getMesh(),
                     getFieldLayout(),
                     GuardCellSizes<Dim>(1),
                     bc_m);
    eg_m.initialize(getMesh(),
                    getFieldLayout(),
                    GuardCellSizes<Dim>(1),
                    vbc_m);

    update();
899 900

//    setGridIsFixed();
gsell's avatar
gsell committed
901 902 903 904 905 906 907 908
}

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

    if(fs_m->hasValidSolver()) {
909 910 911
        //mesh the whole domain
        if(fs_m->getFieldSolverType() == "SAAMG")
            resizeMesh();
gsell's avatar
gsell committed
912 913 914 915 916 917 918 919 920 921 922

        //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;
        gammaz = sqrt(gammaz + 1.0);
923 924
        double scaleFactor = 1;
        // double scaleFactor = Physics::c * getdT();
gsell's avatar
gsell committed
925 926 927 928 929 930 931 932 933
        //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;

934 935 936 937 938
#ifdef DBG_SCALARFIELD
        INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
        ofstream fstr1;
        fstr1.precision(9);

939
        std::string SfileName = OpalData::getInstance()->getInputBasename();
940

941
        std::string rho_fn = std::string("data/") + SfileName + std::string("-rho_scalar-") + std::to_string(fieldDBGStep_m);
942 943 944 945 946 947 948 949 950 951 952 953
        fstr1.open(rho_fn.c_str(), ios::out);
        NDIndex<3> myidx1 = getFieldLayout().getLocalNDIndex();
        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++) {
                    fstr1 << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  rho_m[x][y][z].get() << endl;
                }
            }
        }
        fstr1.close();
        INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif
954

gsell's avatar
gsell committed
955
        // charge density is in rho_m
956
        IpplTimings::startTimer(compPotenTimer_m);
957

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

960
        IpplTimings::stopTimer(compPotenTimer_m);
gsell's avatar
gsell committed
961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980

        //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);

        ofstream fstr2;
        fstr2.precision(9);

981
        std::string phi_fn = std::string("data/") + SfileName + std::string("-phi_scalar-") + std::to_string(fieldDBGStep_m);
982
        fstr2.open(phi_fn.c_str(), ios::out);
gsell's avatar
gsell committed
983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002
        NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
        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++) {
                    fstr2 << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  rho_m[x][y][z].get() << endl;
                }
            }
        }
        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
1003
#ifdef FIELDSTDOUT
1004 1005
        // Immediate debug output:
        // Output potential and e-field along the x-, y-, and z-axes
1006 1007 1008 1009 1010 1011
        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;
1012

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

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

1019
        for (int i=0; i<mz; i++ )
1020
            *gmsg << "Field along z axis Ez = " << eg_m[mx2][my2][i] << " Pot = " << rho_m[mx2][my2][i]  << endl;
1021
#endif
1022

1023
#ifdef DBG_SCALARFIELD
gsell's avatar
gsell committed
1024 1025 1026 1027 1028
        INFOMSG("*** START DUMPING E FIELD ***" << endl);
        //ostringstream oss;
        //MPI_File file;
        //MPI_Status status;
        //MPI_Info fileinfo;
1029
        //MPI_File_open(Ippl::getComm(), "rho_scalar", MPI_MODE_WRONLY | MPI_MODE_CREATE, fileinfo, &file);
gsell's avatar
gsell committed
1030 1031 1032
        ofstream fstr;
        fstr.precision(9);

1033
        std::string e_field = std::string("data/") + SfileName + std::string("-e_field-") + std::to_string(fieldDBGStep_m);
gsell's avatar
gsell committed
1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073
        fstr.open(e_field.c_str(), ios::out);
        NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
        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++) {
                    fstr << x + 1 << " " << y + 1 << " " << z + 1 << " " <<  eg_m[x][y][z].get() << endl;
                }
            }
        }

        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]
         *
         */
        double betaC = sqrt(gammaz * gammaz - 1.0) / gammaz / Physics::c;

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

1074 1075 1076 1077 1078 1079 1080 1081 1082
/**
 * \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()
1083
 * -) Added meanR and quaternion to be handed to the function so that SAAMG solver knows how to rotate the boundary geometry correctly.
1084 1085 1086
 * -) 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?)!
 */
1087
void PartBunch::computeSelfFields_cycl(double gamma) {
1088

1089 1090 1091 1092 1093 1094 1095 1096
    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
1097
    if(fs_m->hasValidSolver()) {
1098 1099 1100
        /// mesh the whole domain
        if(fs_m->getFieldSolverType() == "SAAMG")
            resizeMesh();
1101

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

        /// Lorentz transformation