ParallelTTracker.cpp 52.6 KB
Newer Older
gsell's avatar
gsell committed
1
//
2 3
// Class ParallelTTracker
//   OPAL-T tracker.
gsell's avatar
gsell committed
4 5 6
//   The visitor class for tracking particles with time as independent
//   variable.
//
frey_m's avatar
frey_m committed
7
// Copyright (c) 200x - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
8
// All rights reserved
gsell's avatar
gsell committed
9
//
10 11 12 13 14 15 16 17 18
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
gsell's avatar
gsell committed
19
//
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
#include "AbsBeamline/Monitor.h"
52

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

ParallelTTracker::ParallelTTracker(const Beamline &beamline,
                                   const PartData &reference,
                                   bool revBeam,
58 59 60
                                   bool revTrack):
    Tracker(beamline, reference, revBeam, revTrack),
    itsDataSink_m(NULL),
61
    itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getInitialDirection()),
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
    globalEOL_m(false),
    wakeStatus_m(false),
    wakeFunction_m(NULL),
    pathLength_m(0.0),
    zstart_m(0.0),
    dtCurrentTrack_m(0.0),
    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")),
77
    WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
78
    particleMatterStatus_m(false)
snuverink_j's avatar
snuverink_j committed
79
{}
gsell's avatar
gsell committed
80 81

ParallelTTracker::ParallelTTracker(const Beamline &beamline,
frey_m's avatar
frey_m committed
82
                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
83 84 85 86
                                   DataSink &ds,
                                   const PartData &reference,
                                   bool revBeam,
                                   bool revTrack,
87
                                   const std::vector<unsigned long long> &maxSteps,
88
                                   double zstart,
89
                                   const std::vector<double> &zstop,
90
                                   const std::vector<double> &dt):
frey_m's avatar
frey_m committed
91
    Tracker(beamline, bunch, reference, revBeam, revTrack),
92
    itsDataSink_m(&ds),
93
    itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getInitialDirection()),
94 95 96 97 98 99 100 101 102
    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),
103
    emissionSteps_m(std::numeric_limits<unsigned int>::max()),
104 105 106 107 108
    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")),
109
    WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
110
    particleMatterStatus_m(false)
111
{
112 113
    for (unsigned int i = 0; i < zstop.size(); ++ i) {
        stepSizes_m.push_back(dt[i], zstop[i], maxSteps[i]);
114
    }
115

116 117
    stepSizes_m.sortAscendingZStop();
    stepSizes_m.resetIterator();
gsell's avatar
gsell committed
118 119
}

snuverink_j's avatar
snuverink_j committed
120
ParallelTTracker::~ParallelTTracker() {}
gsell's avatar
gsell committed
121

122 123 124
void ParallelTTracker::visitBeamline(const Beamline &bl) {
    const FlaggedBeamline* fbl = static_cast<const FlaggedBeamline*>(&bl);
    if (fbl->getRelativeFlag()) {
125
        OpalBeamline stash(fbl->getOrigin3D(), fbl->getInitialDirection());
126 127 128 129 130 131 132 133 134
        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);
    }
135
}
136

137
void ParallelTTracker::updateRFElement(std::string elName, double maxPhase) {
138 139
    FieldList cavities = itsOpalBeamline_m.getElementByType(ElementBase::RFCAVITY);
    FieldList travelingwaves = itsOpalBeamline_m.getElementByType(ElementBase::TRAVELINGWAVE);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
140
    cavities.insert(cavities.end(), travelingwaves.begin(), travelingwaves.end());
141 142

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

145
            RFCavity *element = static_cast<RFCavity *>((*fit).getElement().get());
146

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
147
            element->setPhasem(maxPhase);
148
            element->setAutophaseVeto();
149

150 151 152
            INFOMSG("Restored cavity phase from the h5 file. Name: " << element->getName() << ", phase: " << maxPhase << " rad" << endl);
            return;
        }
gsell's avatar
gsell committed
153 154 155
    }
}

156 157
void ParallelTTracker::saveCavityPhases() {
    itsDataSink_m->storeCavityInformation();
gsell's avatar
gsell committed
158 159
}

160 161
void ParallelTTracker::restoreCavityPhases() {
    typedef std::vector<MaxPhasesT>::iterator iterator_t;
162

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

173
void ParallelTTracker::execute() {
kraus's avatar
kraus committed
174
    Inform msg("ParallelTTracker ", *gmsg);
175
    OpalData::getInstance()->setInPrepState(true);
176

177 178
    BorisPusher pusher(itsReference);
    const double globalTimeShift = itsBunch_m->weHaveEnergyBins()? OpalData::getInstance()->getGlobalPhaseShift(): 0.0;
179
    OpalData::getInstance()->setGlobalPhaseShift(0.0);
180

181 182
    // the time step needs to be positive in the setup
    itsBunch_m->setdT(std::abs(itsBunch_m->getdT()));
183
    dtCurrentTrack_m = itsBunch_m->getdT();
184

185
    evenlyDistributeParticles();
186

187
    if (OpalData::getInstance()->hasPriorTrack() || OpalData::getInstance()->inRestartRun()) {
188
        OpalData::getInstance()->setOpenMode(OpalData::OPENMODE::APPEND);
189
    }
190

191
    prepareSections();
192

193 194
    double minTimeStep = stepSizes_m.getMinTimeStep();
    unsigned long long totalNumSteps = stepSizes_m.getNumStepsFinestResolution();
195

196
    itsOpalBeamline_m.activateElements();
197

198 199
    if ( OpalData::getInstance()->hasPriorTrack() ||
         OpalData::getInstance()->inRestartRun()) {
200

201 202
        pathLength_m = itsBunch_m->get_sPos();
        zstart_m = pathLength_m;
203

204 205
        restoreCavityPhases();
    } else {
206

207 208 209 210 211 212
        double momentum = euclidean_norm(itsBunch_m->get_pmean_Distribution());
        CoordinateSystemTrafo beamlineToLab = itsOpalBeamline_m.getCSTrafoLab2Local().inverted();
        itsBunch_m->toLabTrafo_m = beamlineToLab;

        itsBunch_m->RefPartR_m = beamlineToLab.transformTo(Vector_t(0.0));
        itsBunch_m->RefPartP_m = beamlineToLab.rotateTo(momentum * Vector_t(0, 0, 1));
213

214 215
        if (itsBunch_m->getTotalNum() > 0) {
            if (!itsOpalBeamline_m.containsSource()) {
216 217
                momentum = itsReference.getP() / itsBunch_m->getM();
                itsBunch_m->RefPartP_m = beamlineToLab.rotateTo(momentum * Vector_t(0, 0, 1));
218
            }
219

220 221 222
            if (zstart_m > pathLength_m) {
                findStartPosition(pusher);
            }
223

224
            itsBunch_m->set_sPos(pathLength_m);
225 226 227
        }
    }

228 229 230 231 232 233 234 235 236
    stepSizes_m.advanceToPos(zstart_m);
    if (back_track) {
        itsBunch_m->setdT(-std::abs(itsBunch_m->getdT()));
        stepSizes_m.reverseDirection();
        if (pathLength_m < stepSizes_m.getZStop()) {
            ++ stepSizes_m;
        }
    }

237 238 239 240
    Vector_t rmin(0.0), rmax(0.0);
    if (itsBunch_m->getTotalNum() > 0) {
        itsBunch_m->get_bounds(rmin, rmax);
    }
241

242
    OrbitThreader oth(itsReference,
243 244
                      itsBunch_m->RefPartR_m,
                      itsBunch_m->RefPartP_m,
245 246 247
                      pathLength_m,
                      -rmin(2),
                      itsBunch_m->getT(),
248
                      (back_track? -minTimeStep: minTimeStep),
249
                      totalNumSteps,
250
                      stepSizes_m.getFinalZStop() + 2 * rmax(2),
251
                      itsOpalBeamline_m);
252

253
    oth.execute();
254

255
    saveCavityPhases();
256

257
    numParticlesInSimulation_m = itsBunch_m->getTotalNum();
258

259
    setTime();
260

261 262
    double time = itsBunch_m->getT() - globalTimeShift;
    itsBunch_m->setT(time);
263

264
    *gmsg << level1 << *itsBunch_m << endl;
265

266 267
    unsigned long long step = itsBunch_m->getGlobalTrackStep();
    OPALTimer::Timer myt1;
268
    *gmsg << "Track start at: " << myt1.time() << ", t= " << Util::getTimeString(time) << "; "
269 270
          << "zstart at: " << Util::getLengthString(pathLength_m)
          << endl;
271

272
    prepareEmission();
273

274 275
    *gmsg << level1
          << "Executing ParallelTTracker, initial dt= " << Util::getTimeString(itsBunch_m->getdT()) << ";\n"
276
          << "max integration steps " << stepSizes_m.getMaxSteps() << ", next step= " << step << endl;
277

278
    setOptionalVariables();
279

280 281 282 283
    globalEOL_m = false;
    wakeStatus_m = false;
    deletedParticles_m = false;
    OpalData::getInstance()->setInPrepState(false);
284 285 286 287 288 289 290 291 292 293
    while (!stepSizes_m.reachedEnd()) {
        unsigned long long trackSteps = stepSizes_m.getNumSteps() + step;
        dtCurrentTrack_m = stepSizes_m.getdT();
        changeDT(back_track);

        for (; step < trackSteps; ++ step) {
            Vector_t rmin(0.0), rmax(0.0);
            if (itsBunch_m->getTotalNum() > 0) {
                itsBunch_m->get_bounds(rmin, rmax);
            }
294

295
            timeIntegration1(pusher);
296

297 298
            itsBunch_m->Ef = Vector_t(0.0);
            itsBunch_m->Bf = Vector_t(0.0);
299

300
            computeSpaceChargeFields(step);
301

302
            selectDT(back_track);
303
            emitParticles(step);
304
            selectDT(back_track);
305

306
            computeExternalFields(oth);
307

308
            timeIntegration2(pusher);
309

310
            itsBunch_m->incrementT();
311

312 313
            if (itsBunch_m->getT() > 0.0 ||
                itsBunch_m->getdT() < 0.0) {
314
                updateReference(pusher);
315
            }
316

317 318 319
            if (deletedParticles_m) {
                evenlyDistributeParticles();
                deletedParticles_m = false;
320
            }
321
            itsBunch_m->set_sPos(pathLength_m);
322

323
            if (hasEndOfLineReached()) break;
324

325 326 327
            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);
328

329
            itsBunch_m->incTrackSteps();
330

331
            double beta = euclidean_norm(itsBunch_m->RefPartP_m / Util::getGamma(itsBunch_m->RefPartP_m));
332 333 334 335
            double driftPerTimeStep = std::abs(itsBunch_m->getdT()) * Physics::c * beta;
            if (std::abs(stepSizes_m.getZStop() - pathLength_m) < 0.5 * driftPerTimeStep) {
                break;
            }
336
        }
337

338 339
        if (globalEOL_m)
            break;
340

341
        ++ stepSizes_m;
342
    }
343

344
    itsBunch_m->set_sPos(pathLength_m);
345

346 347 348
    if (numParticlesInSimulation_m > minBinEmitted_m) {
        numParticlesInSimulation_m = itsBunch_m->getTotalNum();
    }
349

350 351 352
    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);
353

354
    msg << level2 << "Dump phase space of last step" << endl;
355

356
    itsOpalBeamline_m.switchElementsOff();
357

358 359
    OPALTimer::Timer myt3;
    *gmsg << "done executing ParallelTTracker at " << myt3.time() << endl;
360

361
    Monitor::writeStatistics();
gsell's avatar
gsell committed
362

363
    OpalData::getInstance()->setPriorTrack();
364 365 366 367 368 369
}

void ParallelTTracker::prepareSections() {

    itsBeamline_m.accept(*this);
    itsOpalBeamline_m.prepareSections();
370 371 372 373

    itsOpalBeamline_m.compute3DLattice();
    itsOpalBeamline_m.save3DLattice();
    itsOpalBeamline_m.save3DInput();
374 375
}

376
void ParallelTTracker::timeIntegration1(BorisPusher & pusher) {
377

378 379 380
    IpplTimings::startTimer(timeIntegrationTimer1_m);
    pushParticles(pusher);
    IpplTimings::stopTimer(timeIntegrationTimer1_m);
381 382
}

383

384
void ParallelTTracker::timeIntegration2(BorisPusher & pusher) {
385
    /*
386 387 388 389 390 391 392 393 394 395 396 397 398 399 400
      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.

401
    */
402

403
    IpplTimings::startTimer(timeIntegrationTimer2_m);
404 405 406
    kickParticles(pusher);
    //switchElements();
    pushParticles(pusher);
407

408 409 410
    const unsigned int localNum = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum; ++ i) {
        itsBunch_m->dt[i] = itsBunch_m->getdT();
411
    }
412

413 414 415
    IpplTimings::stopTimer(timeIntegrationTimer2_m);
}

416
void ParallelTTracker::selectDT(bool backTrack) {
417

418 419 420
    if (itsBunch_m->getIfBeamEmitting()) {
        double dt = itsBunch_m->getEmissionDeltaT();
        itsBunch_m->setdT(dt);
421 422
    } else {
        double dt = dtCurrentTrack_m;
423
        itsBunch_m->setdT(dt);
gsell's avatar
gsell committed
424
    }
425 426 427
    if (backTrack) {
        itsBunch_m->setdT(-std::abs(itsBunch_m->getdT()));
    }
428
}
gsell's avatar
gsell committed
429

430 431
void ParallelTTracker::changeDT(bool backTrack) {
    selectDT(backTrack);
432 433 434
    const unsigned int localNum = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum; ++ i) {
        itsBunch_m->dt[i] = itsBunch_m->getdT();
435 436
    }
}
437 438

void ParallelTTracker::emitParticles(long long step) {
439 440 441 442
    if (!itsBunch_m->weHaveEnergyBins()) return;

    if (itsBunch_m->getIfBeamEmitting()) {
        CoordinateSystemTrafo gunToLab = itsOpalBeamline_m.getCSTrafoLab2Local().inverted();
443
        CoordinateSystemTrafo refToGun = itsOpalBeamline_m.getCSTrafoLab2Local() * itsBunch_m->toLabTrafo_m;
444

445
        transformBunch(refToGun);
446

447
        itsBunch_m->switchToUnitlessPositions(true);
448 449 450

        Vector_t externalE = Vector_t(0.0);
        Vector_t externalB = Vector_t(0.0);
451 452 453
        itsOpalBeamline_m.getFieldAt(gunToLab.transformTo(Vector_t(0.0)),
                                     gunToLab.rotateTo(Vector_t(0.0)),
                                     itsBunch_m->getT() + 0.5 * itsBunch_m->getdT(),
454 455
                                     externalE,
                                     externalB);
456
        itsBunch_m->emitParticles(externalE(2));
457

458 459 460
        itsBunch_m->updateNumTotal();
        numParticlesInSimulation_m = itsBunch_m->getTotalNum();
        itsBunch_m->switchOffUnitlessPositions(true);
461

462 463
        transformBunch(refToGun.inverted());
    }
464

465 466 467
    if (step > minStepforReBin_m) {
        itsBunch_m->Rebin();
        itsBunch_m->resetInterpolationCache(true);
gsell's avatar
gsell committed
468
    }
469
}
470

adelmann's avatar
adelmann committed
471

472 473
void ParallelTTracker::computeSpaceChargeFields(unsigned long long step) {
    if (numParticlesInSimulation_m <= minBinEmitted_m) return;
474

475
    if (!itsBunch_m->hasFieldSolver()) return;
gsell's avatar
gsell committed
476

477 478 479 480
    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
481

482 483 484 485 486 487
    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();
488

489 490 491
    if (step % repartFreq_m + 1 == repartFreq_m) {
        doBinaryRepartition();
    }
492

493 494 495 496 497 498
    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) {
499

500 501 502 503
            itsBunch_m->setBinCharge(binNumber);
            itsBunch_m->setGlobalMeanR(itsBunch_m->get_centroid());
            itsBunch_m->computeSelfFields(binNumber);
            itsBunch_m->Q = Q_back;
504

505
        }
506

507
    } else {
508 509
        itsBunch_m->setGlobalMeanR(itsBunch_m->get_centroid());
        itsBunch_m->computeSelfFields();
gsell's avatar
gsell committed
510
    }
511

512 513 514 515 516 517
    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]);
    }
518
}
gsell's avatar
gsell committed
519 520


521 522
void ParallelTTracker::computeExternalFields(OrbitThreader &oth) {
    IpplTimings::startTimer(fieldEvaluationTimer_m);
kraus's avatar
kraus committed
523
    Inform msg("ParallelTTracker ", *gmsg);
524
    const unsigned int localNum = itsBunch_m->getLocalNum();
525
    bool locPartOutOfBounds = false, globPartOutOfBounds = false;
526 527 528
    Vector_t rmin(0.0), rmax(0.0);
    if (itsBunch_m->getTotalNum() > 0)
        itsBunch_m->get_bounds(rmin, rmax);
529
    IndexMap::value_t elements;
530

531 532
    try {
        elements = oth.query(pathLength_m + 0.5 * (rmax(2) + rmin(2)), rmax(2) - rmin(2));
snuverink_j's avatar
snuverink_j committed
533
    } catch(IndexMap::OutOfBounds &e) {
534 535 536
        globalEOL_m = true;
        return;
    }
537

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

541 542
    for (; it != end; ++ it) {
        CoordinateSystemTrafo refToLocalCSTrafo = (itsOpalBeamline_m.getMisalignment((*it)) *
543
                                                   (itsOpalBeamline_m.getCSTrafoLab2Local((*it)) * itsBunch_m->toLabTrafo_m));
544
        CoordinateSystemTrafo localToRefCSTrafo = refToLocalCSTrafo.inverted();
545

546
        (*it)->setCurrentSCoordinate(pathLength_m + rmin(2));
547

548 549
        for (unsigned int i = 0; i < localNum; ++ i) {
            if (itsBunch_m->Bin[i] < 0) continue;
550

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

555
            Vector_t localE(0.0), localB(0.0);
556

557 558 559 560
            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;
561 562
                locPartOutOfBounds = true;

563
                continue;
564
            }
565

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

    IpplTimings::stopTimer(fieldEvaluationTimer_m);
574

575 576 577
    computeWakefield(elements);
    computeParticleMatterInteraction(elements, oth);

578 579
    reduce(locPartOutOfBounds, globPartOutOfBounds, OpOrAssign());

580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597
    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();
        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, "
598
            << "remaining " << numParticlesInSimulation_m << " particles" << endl;
599 600 601 602 603 604 605 606 607 608 609 610 611 612 613
    }
}

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) {
614 615
        if ((*it)->hasWake() && !hasWake) {
            IpplTimings::startTimer(WakeFieldTimer_m);
616

617
            hasWake = true;
618

619 620 621 622 623 624 625 626 627 628 629
            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();
630
            }
631

632
            if (!wfInstance) {
633
                throw OpalException("ParallelTTracker::computeWakefield",
634
                                    "empty wake function");
635
            }
636 637 638 639 640

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

643
            if (!itsBunch_m->hasFieldSolver()) itsBunch_m->calcBeamParameters();
644

645 646 647
            Quaternion alignment = getQuaternion(itsBunch_m->get_pmean(), Vector_t(0, 0, 1));
            CoordinateSystemTrafo referenceToBeamCSTrafo(Vector_t(0.0), alignment);
            CoordinateSystemTrafo beamToReferenceCSTrafo = referenceToBeamCSTrafo.inverted();
648

649 650
            for (unsigned int i = 0; i < localNum; ++ i) {
                itsBunch_m->R[i] = referenceToBeamCSTrafo.transformTo(itsBunch_m->R[i]);
651
                itsBunch_m->P[i] = referenceToBeamCSTrafo.rotateTo(itsBunch_m->P[i]);
652
                itsBunch_m->Ef[i] = referenceToBeamCSTrafo.rotateTo(itsBunch_m->Ef[i]);
653
            }
654

frey_m's avatar
frey_m committed
655
            wfInstance->apply(itsBunch_m);
656

657 658
            for (unsigned int i = 0; i < localNum; ++ i) {
                itsBunch_m->R[i] = beamToReferenceCSTrafo.transformTo(itsBunch_m->R[i]);
659
                itsBunch_m->P[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->P[i]);
660
                itsBunch_m->Ef[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->Ef[i]);
661
            }
662 663

            IpplTimings::stopTimer(WakeFieldTimer_m);
664
        }
665
    }
666

667 668 669 670
    if (wakeStatus_m && !hasWake) {
        msg << level2 << "=============== END WAKE CALCULATION ==============" << endl;
        wakeStatus_m = false;
    }
671
}
gsell's avatar
gsell committed
672

673 674
void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elements, OrbitThreader &oth) {
    Inform msg("ParallelTTracker ", *gmsg);
675
    std::set<IndexMap::value_t::value_type> elementsWithParticleMatterInteraction;
snuverink_j's avatar
snuverink_j committed
676
    std::set<ParticleMatterInteractionHandler*> particleMatterinteractionHandlers;
677 678 679 680
    std::pair<double, double> currentRange(0.0, 0.0);

    while (elements.size() > 0) {
        auto it = elements.begin();
681 682
        if ((*it)->hasParticleMatterInteraction()) {
            elementsWithParticleMatterInteraction.insert(*it);
snuverink_j's avatar
snuverink_j committed
683
            particleMatterinteractionHandlers.insert((*it)->getParticleMatterInteraction());
684 685 686 687 688 689 690

            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());
691
        }
692 693

        elements.erase(it);
694
    }
695

696 697 698 699
    if (elementsWithParticleMatterInteraction.size() > 0) {
        std::set<ParticleMatterInteractionHandler*> oldSPHandlers;
        std::vector<ParticleMatterInteractionHandler*> leftBehindSPHandlers, newSPHandlers;
        for (auto it: activeParticleMatterInteractionHandlers_m) {
700 701 702 703
            oldSPHandlers.insert(it);
        }

        leftBehindSPHandlers.resize(std::max(oldSPHandlers.size(),
snuverink_j's avatar
snuverink_j committed
704
                                             particleMatterinteractionHandlers.size()));
705
        auto last = std::set_difference(oldSPHandlers.begin(), oldSPHandlers.end(),
snuverink_j's avatar
snuverink_j committed
706
                                        particleMatterinteractionHandlers.begin(), particleMatterinteractionHandlers.end(),
707 708 709 710 711
                                        leftBehindSPHandlers.begin());
        leftBehindSPHandlers.resize(last - leftBehindSPHandlers.begin());

        for (auto it: leftBehindSPHandlers) {
            if (!it->stillActive()) {
712
                activeParticleMatterInteractionHandlers_m.erase(it);
713 714 715 716
            }
        }

        newSPHandlers.resize(std::max(oldSPHandlers.size(),
717
                                      elementsWithParticleMatterInteraction.size()));
snuverink_j's avatar
snuverink_j committed
718
        last = std::set_difference(particleMatterinteractionHandlers.begin(), particleMatterinteractionHandlers.end(),
719 720 721 722 723
                                   oldSPHandlers.begin(), oldSPHandlers.end(),
                                   newSPHandlers.begin());
        newSPHandlers.resize(last - newSPHandlers.begin());

        for (auto it: newSPHandlers) {
724
            activeParticleMatterInteractionHandlers_m.insert(it);
725 726
        }

727 728 729
        if(!particleMatterStatus_m) {
            msg << level2 << "============== START PARTICLE MATTER INTERACTION CALCULATION =============" << endl;
            particleMatterStatus_m = true;
730
        }
731
    }
732

733
    if (particleMatterStatus_m) {
734 735 736
        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
737
            ParticleMatterInteractionHandler* onlyDegraderWithParticles = NULL;
738
            int degradersWithParticlesCount = 0;
739
            for (auto it: activeParticleMatterInteractionHandlers_m) {
740 741 742 743 744 745 746 747 748
                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
749 750 751 752
            unsigned int localNum = itsBunch_m->getLocalNum();
            unsigned int totalNum = 0;
            reduce(localNum, totalNum, OpAddAssign());
            bool allParticlesInMat = (totalNum == 0 &&
753 754 755 756 757 758 759
                                      degradersWithParticlesCount == 1);

            if (allParticlesInMat) {
                onlyDegraderWithParticles->setFlagAllParticlesIn(true);
                msg << "All particles in degrader" << endl;
            }

760
            auto boundingSphere = itsBunch_m->getLocalBoundingSphere();
761 762
            unsigned int rediffusedParticles = 0;
            unsigned int numEnteredParticles = 0;
763
            for (auto it: activeParticleMatterInteractionHandlers_m) {
764 765
                ElementBase* element = it->getElement();
                CoordinateSystemTrafo refToLocalCSTrafo = (element->getMisalignment() *
766
                                                           (element->getCSTrafoGlobal2Local() * itsBunch_m->toLabTrafo_m));
767 768 769 770 771 772 773
                CoordinateSystemTrafo localToRefCSTrafo