ParallelTTracker.cpp 55.1 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 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
                                   bool revTrack):
    Tracker(beamline, reference, revBeam, revTrack),
    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")),
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 91 92

    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
                                   beamline.getCoordTransformationTo().conjugate());
    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 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
    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")),
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 137 138 139
    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
                                   beamline.getCoordTransformationTo());
    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
    itsOpalBeamline_m.compute3DLattice();
    itsOpalBeamline_m.save3DLattice();
    itsOpalBeamline_m.save3DInput();
236

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

248 249
        numSteps.pop();
        timeStepSizes.pop();
250 251
    }

252
    itsOpalBeamline_m.activateElements();
253

254 255
    if ( OpalData::getInstance()->hasPriorTrack() ||
         OpalData::getInstance()->inRestartRun()) {
256

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

261 262
        pathLength_m = itsBunch_m->get_sPos();
        zstart_m = pathLength_m;
263

264 265 266
        restoreCavityPhases();
    } else {
        RefPartR_m = Vector_t(0.0);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
267
        RefPartP_m = euclidean_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);
268

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

274 275 276
            if (zstart_m > pathLength_m) {
                findStartPosition(pusher);
            }
277

278
            itsBunch_m->set_sPos(pathLength_m);
279 280 281
        }
    }

282 283 284 285 286 287 288 289 290 291 292 293
    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);
294

295
    oth.execute();
296

297
    saveCavityPhases();
298

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

302
    setTime();
303

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

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

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

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

318
    prepareEmission();
319

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

324
    setOptionalVariables();
325

326
    itsBunch_m->toLabTrafo_m = referenceToLabCSTrafo_m;
327

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

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

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

345
            timeIntegration1(pusher);
346

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

350
            computeSpaceChargeFields(step);
351

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

356
            computeExternalFields(oth);
357

358
            timeIntegration2(pusher);
359

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

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

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

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

377
            if (hasEndOfLineReached()) break;
378

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

383
            itsBunch_m->incTrackSteps();
384

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

390 391
        if (globalEOL_m)
            break;
392

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

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

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

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

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

413
    itsOpalBeamline_m.switchElementsOff();
414

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

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

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

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

void ParallelTTracker::prepareSections() {

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

434
void ParallelTTracker::timeIntegration1(BorisPusher & pusher) {
435

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

446
    IpplTimings::stopTimer(timeIntegrationTimer1_m);
447 448
}

449

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

467
    */
468

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

481 482
    //switchElements();
    pushParticles(pusher);
483
#endif
484

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

490 491 492 493
    IpplTimings::stopTimer(timeIntegrationTimer2_m);
}

void ParallelTTracker::selectDT() {
494

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

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

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

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

519
        transformBunch(refToGun);
520

521
        itsBunch_m->switchToUnitlessPositions(true);
522 523 524

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

532 533 534
        itsBunch_m->updateNumTotal();
        numParticlesInSimulation_m = itsBunch_m->getTotalNum();
        itsBunch_m->switchOffUnitlessPositions(true);
535

536 537
        transformBunch(refToGun.inverted());
    }
538

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

adelmann's avatar
adelmann committed
545

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

549
    if (!itsBunch_m->hasFieldSolver()) return;
gsell's avatar
gsell committed
550

551 552 553 554
    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
555

556 557 558 559 560 561
    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();
562

563 564 565
    if (step % repartFreq_m + 1 == repartFreq_m) {
        doBinaryRepartition();
    }
566

567 568 569 570 571 572
    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) {
573

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

579
        }
580

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

586 587 588 589 590 591
    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]);
    }
592
}
gsell's avatar
gsell committed
593 594


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

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

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

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

619
        (*it)->setCurrentSCoordinate(pathLength_m + rmin(2));
620

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

624 625 626
            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];
627

628
            Vector_t localE(0.0), localB(0.0);
629

630 631 632 633
            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;
634 635
                locPartOutOfBounds = true;

636
                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
    computeWakefield(elements);
    computeParticleMatterInteraction(elements, oth);

651 652
    reduce(locPartOutOfBounds, globPartOutOfBounds, OpOrAssign());

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 687
    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) {
688 689
        if ((*it)->hasWake() && !hasWake) {
            IpplTimings::startTimer(WakeFieldTimer_m);
690

691
            hasWake = true;
692

693 694 695 696 697 698 699 700 701 702 703
            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();
704
            }
705

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

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

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

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

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

frey_m's avatar
frey_m committed
729
            wfInstance->apply(itsBunch_m);
730

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

            IpplTimings::stopTimer(WakeFieldTimer_m);
738
        }
739
    }
740

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

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

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

            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());
765
        }
766 767

        elements.erase(it);
768
    }
769

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

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

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

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

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

801 802 803
        if(!particleMaterStatus_m) {
            msg << level2 << "============== START PARTICLE MATER INTERACTION CALCULATION =============" << endl;
            particleMaterStatus_m = true;
804
        }
Christof Metzger-Kraus's avatar