ParallelTTracker.cpp 56.3 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
// ------------------------------------------------------------------------
// $RCSfile: ParallelTTracker.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelTTracker
//   The visitor class for tracking particles with time as independent
//   variable.
//
// ------------------------------------------------------------------------
//
// $Date: 2004/11/12 20:10:11 $
// $Author: adelmann $
//
// ------------------------------------------------------------------------

kraus's avatar
kraus committed
20 21
#include "Algorithms/ParallelTTracker.h"

gsell's avatar
gsell committed
22 23 24 25 26 27
#include <cfloat>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>
28
#include <limits>
adelmann's avatar
adelmann committed
29
#include <cmath>
gsell's avatar
gsell committed
30

31
#include "Algorithms/PartPusher.h"
32 33
#include "Algorithms/OrbitThreader.h"
#include "Algorithms/CavityAutophaser.h"
34 35
#include "Beamlines/Beamline.h"
#include "Lines/Sequence.h"
36

37
#include "Solvers/CSRWakeFunction.hh"
gsell's avatar
gsell committed
38 39 40 41

#include "AbstractObjects/OpalData.h"

#include "BasicActions/Option.h"
42
#include "Utilities/Options.h"
43
#include "Utilities/Util.h"
gsell's avatar
gsell committed
44 45

#include "Distribution/Distribution.h"
46
#include "ValueDefinitions/RealVariable.h"
gsell's avatar
gsell committed
47 48
#include "Utilities/Timer.h"
#include "Utilities/OpalException.h"
49
#include "Solvers/ParticleMatterInteractionHandler.hh"
gsell's avatar
gsell committed
50
#include "Structure/BoundaryGeometry.h"
51 52
#include "Structure/LossDataSink.h"

gsell's avatar
gsell committed
53 54 55 56 57 58 59
class PartData;

using namespace std;

ParallelTTracker::ParallelTTracker(const Beamline &beamline,
                                   const PartData &reference,
                                   bool revBeam,
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
                                   bool revTrack):
    Tracker(beamline, reference, revBeam, revTrack),
    itsBunch_m(NULL),
    itsDataSink_m(NULL),
    itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getCoordTransformationTo()),
    RefPartR_m(0.0),
    RefPartP_m(0.0),
    globalEOL_m(false),
    wakeStatus_m(false),
    wakeFunction_m(NULL),
    pathLength_m(0.0),
    zstart_m(0.0),
    zStop_m(),
    dtCurrentTrack_m(0.0),
    dtAllTracks_m(),
    localTrackSteps_m(),
    minStepforReBin_m(-1),
    minBinEmitted_m(std::numeric_limits<size_t>::max()),
    repartFreq_m(-1),
    emissionSteps_m(std::numeric_limits<unsigned int>::max()),
    numParticlesInSimulation_m(0),
    timeIntegrationTimer1_m(IpplTimings::getTimer("TIntegration1")),
    timeIntegrationTimer2_m(IpplTimings::getTimer("TIntegration2")),
    fieldEvaluationTimer_m(IpplTimings::getTimer("External field eval")),
    BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
85
    WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
86
    particleMaterStatus_m(false),
87
    totalParticlesInSimulation_m(0)
88
    // , logger_m("designPath_" + std::to_string(Ippl::myNode()) + ".dat")
89
{
90 91 92 93

    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
                                   beamline.getCoordTransformationTo().conjugate());
    referenceToLabCSTrafo_m = labToRef.inverted();
94

95
#ifdef OPAL_DKS
96 97
    if (IpplInfo::DKSEnabled)
        setupDKS();
98
#endif
gsell's avatar
gsell committed
99 100 101 102 103 104 105 106
}

ParallelTTracker::ParallelTTracker(const Beamline &beamline,
                                   PartBunch &bunch,
                                   DataSink &ds,
                                   const PartData &reference,
                                   bool revBeam,
                                   bool revTrack,
107
                                   const std::vector<unsigned long long> &maxSteps,
108
                                   double zstart,
109
                                   const std::vector<double> &zstop,
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
                                   const std::vector<double> &dt):
    Tracker(beamline, reference, revBeam, revTrack),
    itsBunch_m(&bunch),
    itsDataSink_m(&ds),
    itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getCoordTransformationTo()),
    RefPartR_m(0.0),
    RefPartP_m(0.0),
    globalEOL_m(false),
    wakeStatus_m(false),
    wakeFunction_m(NULL),
    pathLength_m(0.0),
    zstart_m(zstart),
    dtCurrentTrack_m(0.0),
    minStepforReBin_m(-1),
    minBinEmitted_m(std::numeric_limits<size_t>::max()),
    repartFreq_m(-1),
    emissionSteps_m(numeric_limits<unsigned int>::max()),
    numParticlesInSimulation_m(0),
    timeIntegrationTimer1_m(IpplTimings::getTimer("TIntegration1")),
    timeIntegrationTimer2_m(IpplTimings::getTimer("TIntegration2")),
    fieldEvaluationTimer_m(IpplTimings::getTimer("External field eval")),
    BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
132
    WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
133
    particleMaterStatus_m(false),
134
    totalParticlesInSimulation_m(0)
135
    // , logger_m("designPath_" + std::to_string(Ippl::myNode()) + ".dat")
136
{
137

138 139 140 141
    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
                                   beamline.getCoordTransformationTo());
    referenceToLabCSTrafo_m = labToRef.inverted();

142 143 144 145 146 147 148 149 150
    for (std::vector<unsigned long long>::const_iterator it = maxSteps.begin(); it != maxSteps.end(); ++ it) {
        localTrackSteps_m.push(*it);
    }
    for (std::vector<double>::const_iterator it = dt.begin(); it != dt.end(); ++ it) {
        dtAllTracks_m.push(*it);
    }
    for (std::vector<double>::const_iterator it = zstop.begin(); it != zstop.end(); ++ it) {
        zStop_m.push(*it);
    }
151

152
#ifdef OPAL_DKS
153 154
    if (IpplInfo::DKSEnabled)
        setupDKS();
155
#endif
gsell's avatar
gsell committed
156 157 158
}

ParallelTTracker::~ParallelTTracker() {
159
#ifdef OPAL_DKS
160 161
    if (IpplInfo::DKSEnabled)
        delete dksbase;
162
#endif
gsell's avatar
gsell committed
163 164
}

165 166 167 168 169 170 171 172 173 174 175 176 177
void ParallelTTracker::visitBeamline(const Beamline &bl) {
    const FlaggedBeamline* fbl = static_cast<const FlaggedBeamline*>(&bl);
    if (fbl->getRelativeFlag()) {
        OpalBeamline stash(fbl->getOrigin3D(), fbl->getCoordTransformationTo());
        stash.swap(itsOpalBeamline_m);
        fbl->iterate(*this, false);
        itsOpalBeamline_m.prepareSections();
        itsOpalBeamline_m.compute3DLattice();
        stash.merge(itsOpalBeamline_m);
        stash.swap(itsOpalBeamline_m);
    } else {
        fbl->iterate(*this, false);
    }
178
}
179

180
void ParallelTTracker::updateRFElement(std::string elName, double maxPhase) {
181 182
    FieldList cavities = itsOpalBeamline_m.getElementByType(ElementBase::RFCAVITY);
    FieldList travelingwaves = itsOpalBeamline_m.getElementByType(ElementBase::TRAVELINGWAVE);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
183
    cavities.insert(cavities.end(), travelingwaves.begin(), travelingwaves.end());
184 185

    for (FieldList::iterator fit = cavities.begin(); fit != cavities.end(); ++ fit) {
186 187
        if ((*fit).getElement()->getName() == elName) {

188
            RFCavity *element = static_cast<RFCavity *>((*fit).getElement().get());
189

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
190
            element->setPhasem(maxPhase);
191
            element->setAutophaseVeto();
192

193 194 195
            INFOMSG("Restored cavity phase from the h5 file. Name: " << element->getName() << ", phase: " << maxPhase << " rad" << endl);
            return;
        }
gsell's avatar
gsell committed
196 197 198
    }
}

199 200
void ParallelTTracker::saveCavityPhases() {
    itsDataSink_m->storeCavityInformation();
gsell's avatar
gsell committed
201 202
}

203 204
void ParallelTTracker::restoreCavityPhases() {
    typedef std::vector<MaxPhasesT>::iterator iterator_t;
205

206 207 208 209 210 211
    if (OpalData::getInstance()->hasPriorTrack() ||
        OpalData::getInstance()->inRestartRun()) {
        iterator_t it = OpalData::getInstance()->getFirstMaxPhases();
        iterator_t end = OpalData::getInstance()->getLastMaxPhases();
        for (; it < end; ++ it) {
            updateRFElement((*it).first, (*it).second);
212
        }
213
    }
214 215
}

216
void ParallelTTracker::execute() {
kraus's avatar
kraus committed
217
    Inform msg("ParallelTTracker ", *gmsg);
218

219
    OpalData::getInstance()->setInPrepState(true);
220

221 222
    BorisPusher pusher(itsReference);
    const double globalTimeShift = itsBunch_m->weHaveEnergyBins()? OpalData::getInstance()->getGlobalPhaseShift(): 0.0;
223
    OpalData::getInstance()->setGlobalPhaseShift(0.0);
224

225
    dtCurrentTrack_m = itsBunch_m->getdT();
226

227
    evenlyDistributeParticles();
228

229 230
    if (OpalData::getInstance()->hasPriorTrack() || OpalData::getInstance()->inRestartRun()) {
        Options::openMode = Options::APPEND;
231
    }
232

233
    prepareSections();
234

235 236 237
    itsOpalBeamline_m.compute3DLattice();
    itsOpalBeamline_m.save3DLattice();
    itsOpalBeamline_m.save3DInput();
238

239 240 241 242 243 244 245 246
    std::queue<double> timeStepSizes(dtAllTracks_m);
    std::queue<unsigned long long> numSteps(localTrackSteps_m);
    double minTimeStep = timeStepSizes.front();
    unsigned long long totalNumSteps = 0;
    while (timeStepSizes.size() > 0) {
        if (minTimeStep > timeStepSizes.front()) {
            totalNumSteps = std::ceil(totalNumSteps * minTimeStep / timeStepSizes.front());
            minTimeStep = timeStepSizes.front();
247
        }
248
        totalNumSteps += std::ceil(numSteps.front() * timeStepSizes.front() / minTimeStep);
249

250 251
        numSteps.pop();
        timeStepSizes.pop();
252 253
    }

254
    itsOpalBeamline_m.activateElements();
255

256 257
    if (OpalData::getInstance()->hasPriorTrack() ||
        OpalData::getInstance()->inRestartRun()) {
258

259 260 261
        referenceToLabCSTrafo_m = itsBunch_m->toLabTrafo_m;
        RefPartR_m = referenceToLabCSTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
        RefPartP_m = referenceToLabCSTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
262

263 264
        pathLength_m = itsBunch_m->get_sPos();
        zstart_m = pathLength_m;
265

266 267 268 269
        restoreCavityPhases();
    } else {
        RefPartR_m = Vector_t(0.0);
        RefPartP_m = euclidian_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);
270

271 272 273 274
        if (itsBunch_m->getTotalNum() > 0) {
            if (!itsOpalBeamline_m.containsSource()) {
                RefPartP_m = euclidian_norm(itsBunch_m->get_pmean())  * Vector_t(0, 0, 1);
            }
275

276 277 278
            if (zstart_m > pathLength_m) {
                findStartPosition(pusher);
            }
279

280
            itsBunch_m->set_sPos(pathLength_m);
281 282 283
        }
    }

284 285 286 287 288 289 290 291 292 293 294 295
    Vector_t rmin, rmax;
    itsBunch_m->get_bounds(rmin, rmax);
    OrbitThreader oth(itsReference,
                      referenceToLabCSTrafo_m.transformTo(RefPartR_m),
                      referenceToLabCSTrafo_m.rotateTo(RefPartP_m),
                      pathLength_m,
                      -rmin(2),
                      itsBunch_m->getT(),
                      minTimeStep,
                      totalNumSteps,
                      zStop_m.back() + 2 * rmax(2),
                      itsOpalBeamline_m);
296

297
    oth.execute();
298

299
    saveCavityPhases();
300

301
    numParticlesInSimulation_m = itsBunch_m->getTotalNum();
302
    totalParticlesInSimulation_m = itsBunch_m->getTotalNum();
303

304
    setTime();
305

306 307
    double t = itsBunch_m->getT() - globalTimeShift;
    itsBunch_m->setT(t);
308

309 310
    itsBunch_m->RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
    itsBunch_m->RefPartP_m = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);
311

312
    *gmsg << level1 << *itsBunch_m << endl;
313

314 315 316 317 318
    unsigned long long step = itsBunch_m->getGlobalTrackStep();
    OPALTimer::Timer myt1;
    *gmsg << "Track start at: " << myt1.time() << ", t= " << Util::getTimeString(t) << "; "
          << "zstart at: " << Util::getLengthString(pathLength_m)
          << endl;
319

320
    prepareEmission();
321

322 323 324
    *gmsg << level1
          << "Executing ParallelTTracker, initial dt= " << Util::getTimeString(itsBunch_m->getdT()) << ";\n"
          << "max integration steps " << getMaxSteps(localTrackSteps_m) << ", next step= " << step << endl;
325

326
    setOptionalVariables();
327

328
    itsBunch_m->toLabTrafo_m = referenceToLabCSTrafo_m;
329

330
#ifdef OPAL_DKS
331 332
    if (IpplInfo::DKSEnabled)
        allocateDeviceMemory();
333
#endif
334

335 336 337 338 339 340 341 342 343
    // loggingFrequency_m = floor(1e-11/itsBunch_m->getdT() + 0.5);
    globalEOL_m = false;
    wakeStatus_m = false;
    deletedParticles_m = false;
    OpalData::getInstance()->setInPrepState(false);
    while (localTrackSteps_m.size() > 0) {
        localTrackSteps_m.front() += step;
        dtCurrentTrack_m = dtAllTracks_m.front();
        changeDT();
344

345
        for (; step < localTrackSteps_m.front(); ++step) {
346

347
            timeIntegration1(pusher);
348

349 350
            itsBunch_m->Ef = Vector_t(0.0);
            itsBunch_m->Bf = Vector_t(0.0);
351

352
            computeSpaceChargeFields(step);
353

354 355 356
            selectDT();
            emitParticles(step);
            selectDT();
357

358
            computeExternalFields(oth);
359

360
            timeIntegration2(pusher);
361

362 363
            t += itsBunch_m->getdT();
            itsBunch_m->setT(t);
364

365 366
            if (t > 0.0) {
                updateRefToLabCSTrafo(pusher);
367
            }
368

369 370 371
            if (deletedParticles_m) {
                evenlyDistributeParticles();
                deletedParticles_m = false;
372
            }
373

374 375 376 377
            itsBunch_m->toLabTrafo_m = referenceToLabCSTrafo_m;
            itsBunch_m->RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
            itsBunch_m->RefPartP_m = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);
            itsBunch_m->set_sPos(pathLength_m);
378

379
            if (hasEndOfLineReached()) break;
380

381 382 383
            bool const psDump = ((itsBunch_m->getGlobalTrackStep() % Options::psDumpFreq) + 1 == Options::psDumpFreq);
            bool const statDump = ((itsBunch_m->getGlobalTrackStep() % Options::statDumpFreq) + 1 == Options::statDumpFreq);
            dumpStats(step, psDump, statDump);
384

385
            itsBunch_m->incTrackSteps();
386

387 388
            double driftPerTimeStep = euclidian_norm(itsBunch_m->getdT() * Physics::c * RefPartP_m / Util::getGamma(RefPartP_m));
            if (std::abs(zStop_m.front() - pathLength_m) < 0.5 * driftPerTimeStep)
389 390
                localTrackSteps_m.front() = step;
        }
391

392 393
        if (globalEOL_m)
            break;
394

395 396 397
        dtAllTracks_m.pop();
        localTrackSteps_m.pop();
        zStop_m.pop();
398
    }
399

400 401 402 403
    itsBunch_m->toLabTrafo_m = referenceToLabCSTrafo_m;
    itsBunch_m->RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
    itsBunch_m->RefPartP_m = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);
    itsBunch_m->set_sPos(pathLength_m);
404

405 406 407
    if (numParticlesInSimulation_m > minBinEmitted_m) {
        numParticlesInSimulation_m = itsBunch_m->getTotalNum();
    }
408

409 410 411
    bool const psDump = (((itsBunch_m->getGlobalTrackStep() - 1) % Options::psDumpFreq) + 1 != Options::psDumpFreq);
    bool const statDump = (((itsBunch_m->getGlobalTrackStep() - 1) % Options::statDumpFreq) + 1 != Options::statDumpFreq);
    writePhaseSpace((step + 1), psDump, statDump);
412

413
    msg << level2 << "Dump phase space of last step" << endl;
414

415
    itsOpalBeamline_m.switchElementsOff();
416

417 418
    OPALTimer::Timer myt3;
    *gmsg << "done executing ParallelTTracker at " << myt3.time() << endl;
419

420
#ifdef OPAL_DKS
421 422
    if (IpplInfo::DKSEnabled)
        freeDeviceMemory();
423
#endif
424

425
    LossDataSink::writeStatistics();
gsell's avatar
gsell committed
426

427
    OpalData::getInstance()->setPriorTrack();
428 429 430 431 432 433 434 435
}

void ParallelTTracker::prepareSections() {

    itsBeamline_m.accept(*this);
    itsOpalBeamline_m.prepareSections();
}

436
void ParallelTTracker::timeIntegration1(BorisPusher & pusher) {
437

438
    IpplTimings::startTimer(timeIntegrationTimer1_m);
439 440 441 442 443 444
#ifdef OPAL_DKS
    if (IpplInfo::DKSEnabled)
        pushParticlesDKS();
    else
        pushParticles(pusher);
#else
445
    pushParticles(pusher);
446
#endif
447

448
    IpplTimings::stopTimer(timeIntegrationTimer1_m);
449 450
}

451

452
void ParallelTTracker::timeIntegration2(BorisPusher & pusher) {
453
    /*
454 455 456 457 458 459 460 461 462 463 464 465 466 467 468
      transport and emit particles
      that passed the cathode in the first
      half-step or that would pass it in the
      second half-step.

      to make IPPL and the field solver happy
      make sure that at least 10 particles are emitted

      also remember that node 0 has
      all the particles to be emitted

      this has to be done *after* the calculation of the
      space charges! thereby we neglect space charge effects
      in the very first step of a new-born particle.

469
    */
470

471
    IpplTimings::startTimer(timeIntegrationTimer2_m);
472
#ifdef OPAL_DKS
473 474 475 476 477 478 479
    if (IpplInfo::DKSEnabled) {
        kickParticlesDKS();
        pushParticlesDKS(false);
    } else {
        kickParticles(pusher);
        pushParticles(pusher);
    }
480
#else
481
    kickParticles(pusher);
482

483 484
    //switchElements();
    pushParticles(pusher);
485
#endif
486

487 488 489
    const unsigned int localNum = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum; ++ i) {
        itsBunch_m->dt[i] = itsBunch_m->getdT();
490
    }
491

492 493 494 495
    IpplTimings::stopTimer(timeIntegrationTimer2_m);
}

void ParallelTTracker::selectDT() {
496

497 498 499
    if (itsBunch_m->getIfBeamEmitting()) {
        double dt = itsBunch_m->getEmissionDeltaT();
        itsBunch_m->setdT(dt);
500 501
    } else {
        double dt = dtCurrentTrack_m;
502
        itsBunch_m->setdT(dt);
gsell's avatar
gsell committed
503
    }
504
}
gsell's avatar
gsell committed
505

506 507
void ParallelTTracker::changeDT() {
    selectDT();
508 509 510
    const unsigned int localNum = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum; ++ i) {
        itsBunch_m->dt[i] = itsBunch_m->getdT();
511 512
    }
}
513 514

void ParallelTTracker::emitParticles(long long step) {
515 516 517 518 519
    if (!itsBunch_m->weHaveEnergyBins()) return;

    if (itsBunch_m->getIfBeamEmitting()) {
        CoordinateSystemTrafo gunToLab = itsOpalBeamline_m.getCSTrafoLab2Local().inverted();
        CoordinateSystemTrafo refToGun = itsOpalBeamline_m.getCSTrafoLab2Local() * referenceToLabCSTrafo_m;
520

521
        transformBunch(refToGun);
522

523
        itsBunch_m->switchToUnitlessPositions(true);
524 525 526

        Vector_t externalE = Vector_t(0.0);
        Vector_t externalB = Vector_t(0.0);
527 528 529
        itsOpalBeamline_m.getFieldAt(gunToLab.transformTo(Vector_t(0.0)),
                                     gunToLab.rotateTo(Vector_t(0.0)),
                                     itsBunch_m->getT() + 0.5 * itsBunch_m->getdT(),
530 531
                                     externalE,
                                     externalB);
532
        itsBunch_m->emitParticles(externalE(2));
533

534 535 536
        itsBunch_m->updateNumTotal();
        numParticlesInSimulation_m = itsBunch_m->getTotalNum();
        itsBunch_m->switchOffUnitlessPositions(true);
537

538 539
        transformBunch(refToGun.inverted());
    }
540

541 542 543
    if (step > minStepforReBin_m) {
        itsBunch_m->Rebin();
        itsBunch_m->resetInterpolationCache(true);
gsell's avatar
gsell committed
544
    }
545
}
546

adelmann's avatar
adelmann committed
547

548 549
void ParallelTTracker::computeSpaceChargeFields(unsigned long long step) {
    if (numParticlesInSimulation_m <= minBinEmitted_m) return;
550

551
    if (!itsBunch_m->hasFieldSolver()) return;
gsell's avatar
gsell committed
552

553 554 555 556
    itsBunch_m->calcBeamParameters();
    Quaternion alignment = getQuaternion(itsBunch_m->get_pmean(), Vector_t(0, 0, 1));
    CoordinateSystemTrafo beamToReferenceCSTrafo(Vector_t(0, 0, pathLength_m), alignment.conjugate());
    CoordinateSystemTrafo referenceToBeamCSTrafo = beamToReferenceCSTrafo.inverted();
gsell's avatar
gsell committed
557

558 559 560 561 562 563
    const unsigned int localNum1 = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum1; ++ i) {
        itsBunch_m->R[i] = referenceToBeamCSTrafo.transformTo(itsBunch_m->R[i]);
    }

    itsBunch_m->boundp();
564

565 566 567
    if (step % repartFreq_m + 1 == repartFreq_m) {
        doBinaryRepartition();
    }
568

569 570 571 572 573 574
    if (itsBunch_m->weHaveEnergyBins()) {
        itsBunch_m->calcGammas();
        itsBunch_m->resetInterpolationCache();
        ParticleAttrib<double> Q_back = itsBunch_m->Q;
        for (int binNumber = 0; binNumber < itsBunch_m->getLastEmittedEnergyBin() &&
                 binNumber < itsBunch_m->getNumberOfEnergyBins(); ++binNumber) {
575

576 577 578 579
            itsBunch_m->setBinCharge(binNumber);
            itsBunch_m->setGlobalMeanR(itsBunch_m->get_centroid());
            itsBunch_m->computeSelfFields(binNumber);
            itsBunch_m->Q = Q_back;
580

581
        }
582

583
    } else {
584 585
        itsBunch_m->setGlobalMeanR(itsBunch_m->get_centroid());
        itsBunch_m->computeSelfFields();
gsell's avatar
gsell committed
586
    }
587

588 589 590 591 592 593
    const unsigned int localNum2 = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum2; ++ i) {
        itsBunch_m->R[i] = beamToReferenceCSTrafo.transformTo(itsBunch_m->R[i]);
        itsBunch_m->Ef[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->Ef[i]);
        itsBunch_m->Bf[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->Bf[i]);
    }
594
}
gsell's avatar
gsell committed
595 596


597 598
void ParallelTTracker::computeExternalFields(OrbitThreader &oth) {
    IpplTimings::startTimer(fieldEvaluationTimer_m);
kraus's avatar
kraus committed
599
    Inform msg("ParallelTTracker ", *gmsg);
600
    const unsigned int localNum = itsBunch_m->getLocalNum();
601

602 603 604
    Vector_t rmin, rmax;
    itsBunch_m->get_bounds(rmin, rmax);
    IndexMap::value_t elements;
605

606 607
    try {
        elements = oth.query(pathLength_m + 0.5 * (rmax(2) + rmin(2)), rmax(2) - rmin(2));
snuverink_j's avatar
snuverink_j committed
608
    } catch(IndexMap::OutOfBounds &e) {
609 610 611
        globalEOL_m = true;
        return;
    }
612

613 614
    IndexMap::value_t::const_iterator it = elements.begin();
    const IndexMap::value_t::const_iterator end = elements.end();
615

616 617 618 619
    for (; it != end; ++ it) {
        CoordinateSystemTrafo refToLocalCSTrafo = (itsOpalBeamline_m.getMisalignment((*it)) *
                                                   (itsOpalBeamline_m.getCSTrafoLab2Local((*it)) * referenceToLabCSTrafo_m));
        CoordinateSystemTrafo localToRefCSTrafo = refToLocalCSTrafo.inverted();
620

621
        (*it)->setCurrentSCoordinate(pathLength_m + rmin(2));
622

623 624
        for (unsigned int i = 0; i < localNum; ++ i) {
            if (itsBunch_m->Bin[i] < 0) continue;
625

626 627 628
            itsBunch_m->R[i] = refToLocalCSTrafo.transformTo(itsBunch_m->R[i]);
            itsBunch_m->P[i] = refToLocalCSTrafo.rotateTo(itsBunch_m->P[i]);
            double &dt = itsBunch_m->dt[i];
629

630
            Vector_t localE(0.0), localB(0.0);
631

632 633 634 635 636
            if ((*it)->apply(i, itsBunch_m->getT() + 0.5 * dt, localE, localB)) {
                itsBunch_m->R[i] = localToRefCSTrafo.transformTo(itsBunch_m->R[i]);
                itsBunch_m->P[i] = localToRefCSTrafo.rotateTo(itsBunch_m->P[i]);
                itsBunch_m->Bin[i] = -1;
                continue;
637
            }
638

639 640 641 642
            itsBunch_m->R[i] = localToRefCSTrafo.transformTo(itsBunch_m->R[i]);
            itsBunch_m->P[i] = localToRefCSTrafo.rotateTo(itsBunch_m->P[i]);
            itsBunch_m->Ef[i] += localToRefCSTrafo.rotateTo(localE);
            itsBunch_m->Bf[i] += localToRefCSTrafo.rotateTo(localB);
643
        }
644 645 646
    }

    IpplTimings::stopTimer(fieldEvaluationTimer_m);
647

648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686
    computeWakefield(elements);
    computeParticleMatterInteraction(elements, oth);

    size_t ne = 0;
    bool globPartOutOfBounds = (min(itsBunch_m->Bin) < 0) && (itsBunch_m->getTotalNum() > 1);
    if (globPartOutOfBounds) {
        if (itsBunch_m->hasFieldSolver()) {
            ne = itsBunch_m->boundp_destroyT();
        } else {
            ne = itsBunch_m->destroyT();
        }
        numParticlesInSimulation_m  = itsBunch_m->getTotalNum();
        totalParticlesInSimulation_m -= ne;
        deletedParticles_m = true;
    }

    size_t totalNum = itsBunch_m->getTotalNum();
    if (numParticlesInSimulation_m > minBinEmitted_m || totalNum > minBinEmitted_m) {
        numParticlesInSimulation_m = totalNum;
    }

    if (ne > 0) {
        msg << level1 << "* Deleted " << ne << " particles, "
            << "remaining " << itsBunch_m->getTotalNum() << " particles" << endl;
    }
}

void ParallelTTracker::computeWakefield(IndexMap::value_t &elements) {
    bool hasWake = false;
    WakeFunction *wfInstance;

    Inform msg("ParallelTTracker ", *gmsg);

    const unsigned int localNum = itsBunch_m->getLocalNum();

    IndexMap::value_t::const_iterator it = elements.begin();
    const IndexMap::value_t::const_iterator end = elements.end();

    for (; it != end; ++ it) {
687 688
        if ((*it)->hasWake() && !hasWake) {
            IpplTimings::startTimer(WakeFieldTimer_m);
689

690
            hasWake = true;
691

692 693 694 695 696 697 698 699 700 701 702
            if ((*it)->getWake()->getType() == "CSRWakeFunction" ||
                (*it)->getWake()->getType() == "CSRIGFWakeFunction") {
                if ((*it)->getType() == ElementBase::RBEND ||
                    (*it)->getType() == ElementBase::SBEND) {
                    wfInstance = (*it)->getWake();
                    wakeFunction_m = wfInstance;
                } else {
                    wfInstance = wakeFunction_m;
                }
            } else {
                wfInstance = (*it)->getWake();
703
            }
704

705
            if (!wfInstance) {
706
                throw OpalException("ParallelTTracker::computeWakefield",
707
                                    "empty wake function");
708
            }
709 710 711 712 713

            if (!wakeStatus_m) {
                msg << level2 << "============== START WAKE CALCULATION =============" << endl;
                wfInstance->initialize((*it).get());
                wakeStatus_m = true;
714 715
            }

716
            if (!itsBunch_m->hasFieldSolver()) itsBunch_m->calcBeamParameters();
717

718 719 720
            Quaternion alignment = getQuaternion(itsBunch_m->get_pmean(), Vector_t(0, 0, 1));
            CoordinateSystemTrafo referenceToBeamCSTrafo(Vector_t(0.0), alignment);
            CoordinateSystemTrafo beamToReferenceCSTrafo = referenceToBeamCSTrafo.inverted();
721

722 723 724
            for (unsigned int i = 0; i < localNum; ++ i) {
                itsBunch_m->R[i] = referenceToBeamCSTrafo.transformTo(itsBunch_m->R[i]);
                itsBunch_m->Ef[i] = referenceToBeamCSTrafo.rotateTo(itsBunch_m->Ef[i]);
725
            }
726

727
            wfInstance->apply(*itsBunch_m);
728

729 730 731
            for (unsigned int i = 0; i < localNum; ++ i) {
                itsBunch_m->R[i] = beamToReferenceCSTrafo.transformTo(itsBunch_m->R[i]);
                itsBunch_m->Ef[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->Ef[i]);
732
            }
733 734

            IpplTimings::stopTimer(WakeFieldTimer_m);
735
        }
736
    }
737

738 739 740 741
    if (wakeStatus_m && !hasWake) {
        msg << level2 << "=============== END WAKE CALCULATION ==============" << endl;
        wakeStatus_m = false;
    }
742
}
gsell's avatar
gsell committed
743

744 745 746 747
void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elements, OrbitThreader &oth) {
    Inform msg("ParallelTTracker ", *gmsg);
    IndexMap::value_t::const_iterator it = elements.begin();
    const IndexMap::value_t::const_iterator end = elements.end();
748
    std::set<IndexMap::value_t::value_type> elementsWithParticleMatterInteraction;
snuverink_j's avatar
snuverink_j committed
749
    std::set<ParticleMatterInteractionHandler*> particleMatterinteractionHandlers;
750 751 752 753
    std::pair<double, double> currentRange(0.0, 0.0);

    while (elements.size() > 0) {
        auto it = elements.begin();
754 755
        if ((*it)->hasParticleMatterInteraction()) {
            elementsWithParticleMatterInteraction.insert(*it);
snuverink_j's avatar
snuverink_j committed
756
            particleMatterinteractionHandlers.insert((*it)->getParticleMatterInteraction());
757 758 759 760 761 762 763

            std::pair<double, double> range = oth.getRange(*it, pathLength_m);
            currentRange.first = std::min(currentRange.first, range.first);
            currentRange.second = std::max(currentRange.second, range.second);

            IndexMap::value_t touching = oth.getTouchingElements(range);
            elements.insert(touching.begin(), touching.end());
764
        }
765 766

        elements.erase(it);
767
    }
768

769 770 771 772
    if (elementsWithParticleMatterInteraction.size() > 0) {
        std::set<ParticleMatterInteractionHandler*> oldSPHandlers;
        std::vector<ParticleMatterInteractionHandler*> leftBehindSPHandlers, newSPHandlers;
        for (auto it: activeParticleMatterInteractionHandlers_m) {
773 774 775 776
            oldSPHandlers.insert(it);
        }

        leftBehindSPHandlers.resize(std::max(oldSPHandlers.size(),
snuverink_j's avatar
snuverink_j committed
777
                                             particleMatterinteractionHandlers.size()));
778
        auto last = std::set_difference(oldSPHandlers.begin(), oldSPHandlers.end(),
snuverink_j's avatar
snuverink_j committed
779
                                        particleMatterinteractionHandlers.begin(), particleMatterinteractionHandlers.end(),
780 781 782 783 784
                                        leftBehindSPHandlers.begin());
        leftBehindSPHandlers.resize(last - leftBehindSPHandlers.begin());

        for (auto it: leftBehindSPHandlers) {
            if (!it->stillActive()) {
785
                activeParticleMatterInteractionHandlers_m.erase(it);
786 787 788 789
            }
        }

        newSPHandlers.resize(std::max(oldSPHandlers.size(),
790
                                      elementsWithParticleMatterInteraction.size()));
snuverink_j's avatar
snuverink_j committed
791
        last = std::set_difference(particleMatterinteractionHandlers.begin(), particleMatterinteractionHandlers.end(),
792 793 794 795 796
                                   oldSPHandlers.begin(), oldSPHandlers.end(),
                                   newSPHandlers.begin());
        newSPHandlers.resize(last - newSPHandlers.begin());

        for (auto it: newSPHandlers) {
797
            activeParticleMatterInteractionHandlers_m.insert(it);
798 799
        }

800 801 802
        if(!particleMaterStatus_m) {
            msg << level2 << "============== START PARTICLE MATER INTERACTION CALCULATION =============" << endl;
            particleMaterStatus_m = true;
803
        }
804
    }
805

806
    if (particleMaterStatus_m) {
807 808 809
        do {
            ///all particles in material if max per node is 2 and other degraders have 0 particles
            //check if more than one degrader has particles
810
            ParticleMatterInteractionHandler* onlyDegraderWithParticles = NULL;
811
            int degradersWithParticlesCount = 0;
812
            for (auto it: activeParticleMatterInteractionHandlers_m) {
813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831
                it->setFlagAllParticlesIn(false);
                if (it->getParticlesInMat() > 0) {
                    onlyDegraderWithParticles = it;
                    ++ degradersWithParticlesCount;
                }
            }

            //if max particles per node is 2, and only one degrader has particles set
            //AllParticlesIn for this degrader to true