ParallelTTracker.cpp 55.4 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 32
#include "Algorithms/OrbitThreader.h"
#include "Algorithms/CavityAutophaser.h"
33
#include "Beamlines/Beamline.h"
34
#include "Beamlines/FlaggedBeamline.h"
35
#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
                                   bool revTrack):
    Tracker(beamline, reference, revBeam, revTrack),
    itsDataSink_m(NULL),
63
    itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getInitialDirection()),
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
    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")),
84
    WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
85
    particleMaterStatus_m(false),
86
    totalParticlesInSimulation_m(0)
87
    // , logger_m("designPath_" + std::to_string(Ippl::myNode()) + ".dat")
88
{
89 90

    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
91
                                   beamline.getInitialDirection());
92
    referenceToLabCSTrafo_m = labToRef.inverted();
93

94
#ifdef OPAL_DKS
95 96
    if (IpplInfo::DKSEnabled)
        setupDKS();
97
#endif
gsell's avatar
gsell committed
98 99 100
}

ParallelTTracker::ParallelTTracker(const Beamline &beamline,
frey_m's avatar
frey_m committed
101
                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
102 103 104 105
                                   DataSink &ds,
                                   const PartData &reference,
                                   bool revBeam,
                                   bool revTrack,
106
                                   const std::vector<unsigned long long> &maxSteps,
107
                                   double zstart,
108
                                   const std::vector<double> &zstop,
109
                                   const std::vector<double> &dt):
frey_m's avatar
frey_m committed
110
    Tracker(beamline, bunch, reference, revBeam, revTrack),
111
    itsDataSink_m(&ds),
112
    itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getInitialDirection()),
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
    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")),
130
    WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
131
    particleMaterStatus_m(false),
132
    totalParticlesInSimulation_m(0)
133
    // , logger_m("designPath_" + std::to_string(Ippl::myNode()) + ".dat")
134
{
135

136
    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
137
                                   beamline.getInitialDirection());
138 139
    referenceToLabCSTrafo_m = labToRef.inverted();

140 141 142 143 144 145 146 147 148
    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);
    }
149

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

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

163 164 165 166 167 168 169 170 171 172 173 174 175
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);
    }
176
}
177

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

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

186
            RFCavity *element = static_cast<RFCavity *>((*fit).getElement().get());
187

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

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

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

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

204 205 206 207 208 209
    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);
210
        }
211
    }
212 213
}

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

217
    OpalData::getInstance()->setInPrepState(true);
218

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

223
    dtCurrentTrack_m = itsBunch_m->getdT();
224

225
    evenlyDistributeParticles();
226

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

231
    prepareSections();
232

233 234 235 236 237 238 239 240
    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();
241
        }
242
        totalNumSteps += std::ceil(numSteps.front() * timeStepSizes.front() / minTimeStep);
243

244 245
        numSteps.pop();
        timeStepSizes.pop();
246 247
    }

248
    itsOpalBeamline_m.activateElements();
249

250 251
    if ( OpalData::getInstance()->hasPriorTrack() ||
         OpalData::getInstance()->inRestartRun()) {
252

253 254 255
        referenceToLabCSTrafo_m = itsBunch_m->toLabTrafo_m;
        RefPartR_m = referenceToLabCSTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
        RefPartP_m = referenceToLabCSTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
256

257 258
        pathLength_m = itsBunch_m->get_sPos();
        zstart_m = pathLength_m;
259

260 261
        restoreCavityPhases();
    } else {
262

263
        RefPartR_m = Vector_t(0.0);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
264
        RefPartP_m = euclidean_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);
265

266 267
        if (itsBunch_m->getTotalNum() > 0) {
            if (!itsOpalBeamline_m.containsSource()) {
268
                RefPartP_m = itsReference.getP() / itsBunch_m->getM() * Vector_t(0, 0, 1);
269
            }
270

271 272 273
            if (zstart_m > pathLength_m) {
                findStartPosition(pusher);
            }
274

275
            itsBunch_m->set_sPos(pathLength_m);
276 277 278
        }
    }

279 280 281 282
    Vector_t rmin(0.0), rmax(0.0);
    if (itsBunch_m->getTotalNum() > 0) {
        itsBunch_m->get_bounds(rmin, rmax);
    }
283 284 285 286 287 288 289 290 291 292
    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);
293

294
    oth.execute();
295

296
    saveCavityPhases();
297

298
    numParticlesInSimulation_m = itsBunch_m->getTotalNum();
299
    totalParticlesInSimulation_m = itsBunch_m->getTotalNum();
300

301
    setTime();
302

303 304
    double t = itsBunch_m->getT() - globalTimeShift;
    itsBunch_m->setT(t);
305

306 307
    itsBunch_m->RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
    itsBunch_m->RefPartP_m = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);
308

309
    *gmsg << level1 << *itsBunch_m << endl;
310

311 312 313 314 315
    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;
316

317
    prepareEmission();
318

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

323
    setOptionalVariables();
324

325
    itsBunch_m->toLabTrafo_m = referenceToLabCSTrafo_m;
326

327
#ifdef OPAL_DKS
328 329
    if (IpplInfo::DKSEnabled)
        allocateDeviceMemory();
330
#endif
331

332 333 334 335 336 337 338 339 340
    // 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();
341

342
        for (; step < localTrackSteps_m.front(); ++step) {
343

344
            timeIntegration1(pusher);
345

346 347
            itsBunch_m->Ef = Vector_t(0.0);
            itsBunch_m->Bf = Vector_t(0.0);
348

349
            computeSpaceChargeFields(step);
350

351 352 353
            selectDT();
            emitParticles(step);
            selectDT();
354

355
            computeExternalFields(oth);
356

357
            timeIntegration2(pusher);
358

359 360
            t += itsBunch_m->getdT();
            itsBunch_m->setT(t);
361

362 363
            if (t > 0.0) {
                updateRefToLabCSTrafo(pusher);
364
            }
365

366 367 368
            if (deletedParticles_m) {
                evenlyDistributeParticles();
                deletedParticles_m = false;
369
            }
370

371 372 373 374
            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);
375

376
            if (hasEndOfLineReached()) break;
377

378 379 380
            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);
381

382
            itsBunch_m->incTrackSteps();
383

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
384
            double driftPerTimeStep = euclidean_norm(itsBunch_m->getdT() * Physics::c * RefPartP_m / Util::getGamma(RefPartP_m));
385
            if (std::abs(zStop_m.front() - pathLength_m) < 0.5 * driftPerTimeStep)
386 387
                localTrackSteps_m.front() = step;
        }
388

389 390
        if (globalEOL_m)
            break;
391

392 393 394
        dtAllTracks_m.pop();
        localTrackSteps_m.pop();
        zStop_m.pop();
395
    }
396

397 398 399 400
    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);
401

402 403 404
    if (numParticlesInSimulation_m > minBinEmitted_m) {
        numParticlesInSimulation_m = itsBunch_m->getTotalNum();
    }
405

406 407 408
    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);
409

410
    msg << level2 << "Dump phase space of last step" << endl;
411

412
    itsOpalBeamline_m.switchElementsOff();
413

414 415
    OPALTimer::Timer myt3;
    *gmsg << "done executing ParallelTTracker at " << myt3.time() << endl;
416

417
#ifdef OPAL_DKS
418 419
    if (IpplInfo::DKSEnabled)
        freeDeviceMemory();
420
#endif
421

422
    LossDataSink::writeStatistics();
gsell's avatar
gsell committed
423

424
    OpalData::getInstance()->setPriorTrack();
425 426 427 428 429 430
}

void ParallelTTracker::prepareSections() {

    itsBeamline_m.accept(*this);
    itsOpalBeamline_m.prepareSections();
431 432 433 434

    itsOpalBeamline_m.compute3DLattice();
    itsOpalBeamline_m.save3DLattice();
    itsOpalBeamline_m.save3DInput();
435 436
}

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

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

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

452

453
void ParallelTTracker::timeIntegration2(BorisPusher & pusher) {
454
    /*
455 456 457 458 459 460 461 462 463 464 465 466 467 468 469
      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.

470
    */
471

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

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

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

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

void ParallelTTracker::selectDT() {
497

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

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

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

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

522
        transformBunch(refToGun);
523

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

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

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

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

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

adelmann's avatar
adelmann committed
548

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

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

554 555 556 557
    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
558

559 560 561 562 563 564
    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();
565

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

570 571 572 573 574 575
    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) {
576

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

582
        }
583

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

589 590 591 592 593 594
    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]);
    }
595
}
gsell's avatar
gsell committed
596 597


598 599
void ParallelTTracker::computeExternalFields(OrbitThreader &oth) {
    IpplTimings::startTimer(fieldEvaluationTimer_m);
kraus's avatar
kraus committed
600
    Inform msg("ParallelTTracker ", *gmsg);
601
    const unsigned int localNum = itsBunch_m->getLocalNum();
602
    bool locPartOutOfBounds = false, globPartOutOfBounds = false;
603 604 605
    Vector_t rmin(0.0), rmax(0.0);
    if (itsBunch_m->getTotalNum() > 0)
        itsBunch_m->get_bounds(rmin, rmax);
606
    IndexMap::value_t elements;
607

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

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

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

623
        (*it)->setCurrentSCoordinate(pathLength_m + rmin(2));
624

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

628 629 630
            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];
631

632
            Vector_t localE(0.0), localB(0.0);
633

634 635 636 637
            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;
638 639
                locPartOutOfBounds = true;

640
                continue;
641
            }
642

643 644 645 646
            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);
647
        }
648 649 650
    }

    IpplTimings::stopTimer(fieldEvaluationTimer_m);
651

652 653 654
    computeWakefield(elements);
    computeParticleMatterInteraction(elements, oth);

655 656
    reduce(locPartOutOfBounds, globPartOutOfBounds, OpOrAssign());

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 687 688 689 690 691
    size_t ne = 0;
    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) {
692 693
        if ((*it)->hasWake() && !hasWake) {
            IpplTimings::startTimer(WakeFieldTimer_m);
694

695
            hasWake = true;
696

697 698 699 700 701 702 703 704 705 706 707
            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();
708
            }
709

710
            if (!wfInstance) {
711
                throw OpalException("ParallelTTracker::computeWakefield",
712
                                    "empty wake function");
713
            }
714 715 716 717 718

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

721
            if (!itsBunch_m->hasFieldSolver()) itsBunch_m->calcBeamParameters();
722

723 724 725
            Quaternion alignment = getQuaternion(itsBunch_m->get_pmean(), Vector_t(0, 0, 1));
            CoordinateSystemTrafo referenceToBeamCSTrafo(Vector_t(0.0), alignment);
            CoordinateSystemTrafo beamToReferenceCSTrafo = referenceToBeamCSTrafo.inverted();
726

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

frey_m's avatar
frey_m committed
733
            wfInstance->apply(itsBunch_m);
734

735 736
            for (unsigned int i = 0; i < localNum; ++ i) {
                itsBunch_m->R[i] = beamToReferenceCSTrafo.transformTo(itsBunch_m->R[i]);
737
                itsBunch_m->P[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->P[i]);
738
                itsBunch_m->Ef[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->Ef[i]);
739
            }
740 741

            IpplTimings::stopTimer(WakeFieldTimer_m);
742
        }
743
    }
744

745 746 747 748
    if (wakeStatus_m && !hasWake) {
        msg << level2 << "=============== END WAKE CALCULATION ==============" << endl;
        wakeStatus_m = false;
    }
749
}
gsell's avatar
gsell committed
750

751 752
void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elements, OrbitThreader &oth) {
    Inform msg("ParallelTTracker ", *gmsg);
753
    std::set<IndexMap::value_t::value_type> elementsWithParticleMatterInteraction;
snuverink_j's avatar
snuverink_j committed
754
    std::set<ParticleMatterInteractionHandler*> particleMatterinteractionHandlers;
755 756 757 758
    std::pair<double, double> currentRange(0.0, 0.0);

    while (elements.size() > 0) {
        auto it = elements.begin();
759 760
        if ((*it)->hasParticleMatterInteraction()) {
            elementsWithParticleMatterInteraction.insert(*it);
snuverink_j's avatar
snuverink_j committed
761
            particleMatterinteractionHandlers.insert((*it)->getParticleMatterInteraction());
762 763 764 765 766 767 768

            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());
769
        }
770 771

        elements.erase(it);
772
    }
773

774 775 776 777
    if (elementsWithParticleMatterInteraction.size() > 0) {
        std::set<ParticleMatterInteractionHandler*> oldSPHandlers;
        std::vector<ParticleMatterInteractionHandler*> leftBehindSPHandlers, newSPHandlers;
        for (auto it: activeParticleMatterInteractionHandlers_m) {
778 779 780 781
            oldSPHandlers.insert(it);
        }

        leftBehindSPHandlers.resize(std::max(oldSPHandlers.size(),
snuverink_j's avatar
snuverink_j committed
782
                                             particleMatterinteractionHandlers.size()));
783
        auto last = std::set_difference(oldSPHandlers.begin(), oldSPHandlers.end(),
snuverink_j's avatar
snuverink_j committed
784
                                        particleMatterinteractionHandlers.begin(), particleMatterinteractionHandlers.end(),
785 786 787 788 789
                                        leftBehindSPHandlers.begin());
        leftBehindSPHandlers.resize(last - leftBehindSPHandlers.begin());

        for (auto it: leftBehindSPHandlers) {
            if (!it->stillActive()) {
790
                activeParticleMatterInte