ParallelTTracker.cpp 56.7 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// ------------------------------------------------------------------------
// $RCSfile: ParallelTTracker.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelTTracker
//   The visitor class for tracking particles with time as independent
//   variable.
//
// ------------------------------------------------------------------------
//
// $Date: 2004/11/12 20:10:11 $
// $Author: adelmann $
//
// ------------------------------------------------------------------------

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

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

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

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

#include "AbstractObjects/OpalData.h"

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

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

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

using namespace std;

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

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

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

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

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

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

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

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

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

180
void ParallelTTracker::updateRFElement(std::string elName, double maxPhase) {
181
    /**
182 183 184 185 186 187 188 189
       The maximum phase is added to the nominal phase of
       the element. This is done on all nodes except node 0 where
       the Autophase took place.
    */
    FieldList cavities = itsOpalBeamline_m.getElementByType(ElementBase::RFCAVITY);
    FieldList travelingwaves = itsOpalBeamline_m.getElementByType(ElementBase::TRAVELINGWAVE);

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

192 193
            RFCavity *element = static_cast<RFCavity *>((*fit).getElement().get());
            double phase  =  element->getPhasem();
194

195 196
            element->setPhasem(phase + maxPhase);
            element->setAutophaseVeto();
197

198 199
            INFOMSG("Restored cavity phase from the h5 file. Name: " << element->getName() << ", phase: " << maxPhase << " rad" << endl);
            return;
200
        }
gsell's avatar
gsell committed
201
    }
202 203
    for (FieldList::iterator fit = travelingwaves.begin(); fit != travelingwaves.end(); ++ fit) {
        if ((*fit).getElement()->getName() == elName) {
204

205 206
            TravelingWave *element = static_cast<TravelingWave *>((*fit).getElement().get());
            double phase  =  element->getPhasem();
207

208 209
            element->setPhasem(phase + maxPhase);
            element->setAutophaseVeto();
210

211 212 213
            INFOMSG("Restored cavity phase from the h5 file. Name: " << element->getName() << ", phase: " << maxPhase << " rad" << endl);
            return;
        }
gsell's avatar
gsell committed
214 215 216
    }
}

217 218
void ParallelTTracker::saveCavityPhases() {
    itsDataSink_m->storeCavityInformation();
gsell's avatar
gsell committed
219 220
}

221 222
void ParallelTTracker::restoreCavityPhases() {
    typedef std::vector<MaxPhasesT>::iterator iterator_t;
223

224 225 226 227 228 229
    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);
230
        }
231
    }
232 233
}

234
void ParallelTTracker::execute() {
kraus's avatar
kraus committed
235
    Inform msg("ParallelTTracker ", *gmsg);
236

237
    OpalData::getInstance()->setInPrepState(true);
238

239 240
    BorisPusher pusher(itsReference);
    const double globalTimeShift = itsBunch_m->weHaveEnergyBins()? OpalData::getInstance()->getGlobalPhaseShift(): 0.0;
241
    OpalData::getInstance()->setGlobalPhaseShift(0.0);
242

243
    dtCurrentTrack_m = itsBunch_m->getdT();
244

245
    evenlyDistributeParticles();
246

247 248
    if (OpalData::getInstance()->hasPriorTrack() || OpalData::getInstance()->inRestartRun()) {
        Options::openMode = Options::APPEND;
249
    }
250

251
    prepareSections();
252

253 254 255
    itsOpalBeamline_m.compute3DLattice();
    itsOpalBeamline_m.save3DLattice();
    itsOpalBeamline_m.save3DInput();
256

257 258 259 260 261 262 263 264
    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();
265
        }
266
        totalNumSteps += std::ceil(numSteps.front() * timeStepSizes.front() / minTimeStep);
267

268 269
        numSteps.pop();
        timeStepSizes.pop();
270 271
    }

272
    itsOpalBeamline_m.activateElements();
273

274 275
    if (OpalData::getInstance()->hasPriorTrack() ||
        OpalData::getInstance()->inRestartRun()) {
276

277 278 279
        referenceToLabCSTrafo_m = itsBunch_m->toLabTrafo_m;
        RefPartR_m = referenceToLabCSTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
        RefPartP_m = referenceToLabCSTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
280

281 282
        pathLength_m = itsBunch_m->get_sPos();
        zstart_m = pathLength_m;
283

284 285 286 287
        restoreCavityPhases();
    } else {
        RefPartR_m = Vector_t(0.0);
        RefPartP_m = euclidian_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);
288

289 290 291 292
        if (itsBunch_m->getTotalNum() > 0) {
            if (!itsOpalBeamline_m.containsSource()) {
                RefPartP_m = euclidian_norm(itsBunch_m->get_pmean())  * Vector_t(0, 0, 1);
            }
293

294 295 296
            if (zstart_m > pathLength_m) {
                findStartPosition(pusher);
            }
297

298
            itsBunch_m->set_sPos(pathLength_m);
299 300 301
        }
    }

302 303 304 305 306 307 308 309 310 311 312 313
    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);
314

315
    oth.execute();
316

317
    saveCavityPhases();
318

319
    numParticlesInSimulation_m = itsBunch_m->getTotalNum();
320
    totalParticlesInSimulation_m = itsBunch_m->getTotalNum();
321

322
    setTime();
323

324 325
    double t = itsBunch_m->getT() - globalTimeShift;
    itsBunch_m->setT(t);
326

327 328
    itsBunch_m->RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
    itsBunch_m->RefPartP_m = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);
329

330
    *gmsg << level1 << *itsBunch_m << endl;
331

332 333 334 335 336
    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;
337

338
    prepareEmission();
339

340 341 342
    *gmsg << level1
          << "Executing ParallelTTracker, initial dt= " << Util::getTimeString(itsBunch_m->getdT()) << ";\n"
          << "max integration steps " << getMaxSteps(localTrackSteps_m) << ", next step= " << step << endl;
343

344
    setOptionalVariables();
345

346
    itsBunch_m->toLabTrafo_m = referenceToLabCSTrafo_m;
347

348
#ifdef OPAL_DKS
349 350
    if (IpplInfo::DKSEnabled)
        allocateDeviceMemory();
351
#endif
352

353 354 355 356 357 358 359 360 361
    // 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();
362

363
        for (; step < localTrackSteps_m.front(); ++step) {
364

365
            timeIntegration1(pusher);
366

367 368
            itsBunch_m->Ef = Vector_t(0.0);
            itsBunch_m->Bf = Vector_t(0.0);
369

370
            computeSpaceChargeFields(step);
371

372 373 374
            selectDT();
            emitParticles(step);
            selectDT();
375

376
            computeExternalFields(oth);
377

378
            timeIntegration2(pusher);
379

380 381
            t += itsBunch_m->getdT();
            itsBunch_m->setT(t);
382

383 384
            if (t > 0.0) {
                updateRefToLabCSTrafo(pusher);
385
            }
386

387 388 389
            if (deletedParticles_m) {
                evenlyDistributeParticles();
                deletedParticles_m = false;
390
            }
391 392 393 394
            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);
395

396
            if (hasEndOfLineReached()) break;
397

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

402
            itsBunch_m->incTrackSteps();
403

404 405
            double driftPerTimeStep = euclidian_norm(itsBunch_m->getdT() * Physics::c * RefPartP_m / Util::getGamma(RefPartP_m));
            if (std::abs(zStop_m.front() - pathLength_m) < 0.5 * driftPerTimeStep)
406 407
                localTrackSteps_m.front() = step;
        }
408

409 410
        if (globalEOL_m)
            break;
411

412 413 414
        dtAllTracks_m.pop();
        localTrackSteps_m.pop();
        zStop_m.pop();
415
    }
416

417 418 419 420
    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);
421

422 423 424
    if (numParticlesInSimulation_m > minBinEmitted_m) {
        numParticlesInSimulation_m = itsBunch_m->getTotalNum();
    }
425

426 427 428
    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);
429

430
    msg << level2 << "Dump phase space of last step" << endl;
431

432
    itsOpalBeamline_m.switchElementsOff();
433

434 435
    OPALTimer::Timer myt3;
    *gmsg << "done executing ParallelTTracker at " << myt3.time() << endl;
436

437
#ifdef OPAL_DKS
438 439
    if (IpplInfo::DKSEnabled)
        freeDeviceMemory();
440
#endif
441

442
    LossDataSink::writeStatistics();
gsell's avatar
gsell committed
443

444
    OpalData::getInstance()->setPriorTrack();
445 446 447 448 449 450 451 452
}

void ParallelTTracker::prepareSections() {

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

453
void ParallelTTracker::timeIntegration1(BorisPusher & pusher) {
454

455
    IpplTimings::startTimer(timeIntegrationTimer1_m);
456 457 458 459 460 461
#ifdef OPAL_DKS
    if (IpplInfo::DKSEnabled)
        pushParticlesDKS();
    else
        pushParticles(pusher);
#else
462
    pushParticles(pusher);
463
#endif
464

465
    IpplTimings::stopTimer(timeIntegrationTimer1_m);
466 467
}

468

469
void ParallelTTracker::timeIntegration2(BorisPusher & pusher) {
470
    /*
471 472 473 474 475 476 477 478 479 480 481 482 483 484 485
      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.

486
    */
487

488
    IpplTimings::startTimer(timeIntegrationTimer2_m);
489
#ifdef OPAL_DKS
490 491 492 493 494 495 496
    if (IpplInfo::DKSEnabled) {
        kickParticlesDKS();
        pushParticlesDKS(false);
    } else {
        kickParticles(pusher);
        pushParticles(pusher);
    }
497
#else
498
    kickParticles(pusher);
499

500 501
    //switchElements();
    pushParticles(pusher);
502
#endif
503

504 505 506
    const unsigned int localNum = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum; ++ i) {
        itsBunch_m->dt[i] = itsBunch_m->getdT();
507
    }
508

509 510 511 512
    IpplTimings::stopTimer(timeIntegrationTimer2_m);
}

void ParallelTTracker::selectDT() {
513

514 515 516
    if (itsBunch_m->getIfBeamEmitting()) {
        double dt = itsBunch_m->getEmissionDeltaT();
        itsBunch_m->setdT(dt);
517 518
    } else {
        double dt = dtCurrentTrack_m;
519
        itsBunch_m->setdT(dt);
gsell's avatar
gsell committed
520
    }
521
}
gsell's avatar
gsell committed
522

523 524
void ParallelTTracker::changeDT() {
    selectDT();
525 526 527
    const unsigned int localNum = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum; ++ i) {
        itsBunch_m->dt[i] = itsBunch_m->getdT();
528 529
    }
}
530 531

void ParallelTTracker::emitParticles(long long step) {
532 533 534 535 536
    if (!itsBunch_m->weHaveEnergyBins()) return;

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

538
        transformBunch(refToGun);
539

540
        itsBunch_m->switchToUnitlessPositions(true);
541 542 543

        Vector_t externalE = Vector_t(0.0);
        Vector_t externalB = Vector_t(0.0);
544 545 546
        itsOpalBeamline_m.getFieldAt(gunToLab.transformTo(Vector_t(0.0)),
                                     gunToLab.rotateTo(Vector_t(0.0)),
                                     itsBunch_m->getT() + 0.5 * itsBunch_m->getdT(),
547 548
                                     externalE,
                                     externalB);
549
        itsBunch_m->emitParticles(externalE(2));
550

551 552 553
        itsBunch_m->updateNumTotal();
        numParticlesInSimulation_m = itsBunch_m->getTotalNum();
        itsBunch_m->switchOffUnitlessPositions(true);
554

555 556
        transformBunch(refToGun.inverted());
    }
557

558 559 560
    if (step > minStepforReBin_m) {
        itsBunch_m->Rebin();
        itsBunch_m->resetInterpolationCache(true);
gsell's avatar
gsell committed
561
    }
562
}
563

adelmann's avatar
adelmann committed
564

565 566
void ParallelTTracker::computeSpaceChargeFields(unsigned long long step) {
    if (numParticlesInSimulation_m <= minBinEmitted_m) return;
567

568
    if (!itsBunch_m->hasFieldSolver()) return;
gsell's avatar
gsell committed
569

570 571 572 573
    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
574

575 576 577 578 579 580
    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();
581

582 583 584
    if (step % repartFreq_m + 1 == repartFreq_m) {
        doBinaryRepartition();
    }
585

586 587 588 589 590 591
    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) {
592

593 594 595 596
            itsBunch_m->setBinCharge(binNumber);
            itsBunch_m->setGlobalMeanR(itsBunch_m->get_centroid());
            itsBunch_m->computeSelfFields(binNumber);
            itsBunch_m->Q = Q_back;
597

598
        }
599

600
    } else {
601 602
        itsBunch_m->setGlobalMeanR(itsBunch_m->get_centroid());
        itsBunch_m->computeSelfFields();
gsell's avatar
gsell committed
603
    }
604

605 606 607 608 609 610
    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]);
    }
611
}
gsell's avatar
gsell committed
612 613


614 615
void ParallelTTracker::computeExternalFields(OrbitThreader &oth) {
    IpplTimings::startTimer(fieldEvaluationTimer_m);
kraus's avatar
kraus committed
616
    Inform msg("ParallelTTracker ", *gmsg);
617
    const unsigned int localNum = itsBunch_m->getLocalNum();
618

619 620 621
    Vector_t rmin, rmax;
    itsBunch_m->get_bounds(rmin, rmax);
    IndexMap::value_t elements;
622

623 624
    try {
        elements = oth.query(pathLength_m + 0.5 * (rmax(2) + rmin(2)), rmax(2) - rmin(2));
snuverink_j's avatar
snuverink_j committed
625
    } catch(IndexMap::OutOfBounds &e) {
626 627 628
        globalEOL_m = true;
        return;
    }
629

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

633 634 635 636
    for (; it != end; ++ it) {
        CoordinateSystemTrafo refToLocalCSTrafo = (itsOpalBeamline_m.getMisalignment((*it)) *
                                                   (itsOpalBeamline_m.getCSTrafoLab2Local((*it)) * referenceToLabCSTrafo_m));
        CoordinateSystemTrafo localToRefCSTrafo = refToLocalCSTrafo.inverted();
637

638
        (*it)->setCurrentSCoordinate(pathLength_m + rmin(2));
639

640 641
        for (unsigned int i = 0; i < localNum; ++ i) {
            if (itsBunch_m->Bin[i] < 0) continue;
642

643 644 645
            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];
646

647
            Vector_t localE(0.0), localB(0.0);
648

649 650 651 652 653
            if ((*it)->apply(i, itsBunch_m->getT() + 0.5 * dt, localE, localB)) {
                itsBunch_m->R[i] = localToRefCSTrafo.transformTo(itsBunch_m->R[i]);
                itsBunch_m->P[i] = localToRefCSTrafo.rotateTo(itsBunch_m->P[i]);
                itsBunch_m->Bin[i] = -1;
                continue;
654
            }
655

656 657 658 659
            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);
660
        }
661 662 663
    }

    IpplTimings::stopTimer(fieldEvaluationTimer_m);
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 692 693 694 695 696 697 698 699 700 701 702 703
    computeWakefield(elements);
    computeParticleMatterInteraction(elements, oth);

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

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

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

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

    Inform msg("ParallelTTracker ", *gmsg);

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

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

    for (; it != end; ++ it) {
704 705
        if ((*it)->hasWake() && !hasWake) {
            IpplTimings::startTimer(WakeFieldTimer_m);
706

707
            hasWake = true;
708

709 710 711 712 713 714 715 716 717 718 719
            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();
720
            }
721

722
            if (!wfInstance) {
723
                throw OpalException("ParallelTTracker::computeWakefield",
724
                                    "empty wake function");
725
            }
726 727 728 729 730

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

733
            if (!itsBunch_m->hasFieldSolver()) itsBunch_m->calcBeamParameters();
734

735 736 737
            Quaternion alignment = getQuaternion(itsBunch_m->get_pmean(), Vector_t(0, 0, 1));
            CoordinateSystemTrafo referenceToBeamCSTrafo(Vector_t(0.0), alignment);
            CoordinateSystemTrafo beamToReferenceCSTrafo = referenceToBeamCSTrafo.inverted();
738

739 740 741
            for (unsigned int i = 0; i < localNum; ++ i) {
                itsBunch_m->R[i] = referenceToBeamCSTrafo.transformTo(itsBunch_m->R[i]);
                itsBunch_m->Ef[i] = referenceToBeamCSTrafo.rotateTo(itsBunch_m->Ef[i]);
742
            }
743

744
            wfInstance->apply(*itsBunch_m);
745

746 747 748
            for (unsigned int i = 0; i < localNum; ++ i) {
                itsBunch_m->R[i] = beamToReferenceCSTrafo.transformTo(itsBunch_m->R[i]);
                itsBunch_m->Ef[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->Ef[i]);
749
            }
750 751

            IpplTimings::stopTimer(WakeFieldTimer_m);
752
        }
753
    }
754

755 756 757 758
    if (wakeStatus_m && !hasWake) {
        msg << level2 << "=============== END WAKE CALCULATION ==============" << endl;
        wakeStatus_m = false;
    }
759
}
gsell's avatar
gsell committed
760

761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780
void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elements, OrbitThreader &oth) {
    Inform msg("ParallelTTracker ", *gmsg);
    IndexMap::value_t::const_iterator it = elements.begin();
    const IndexMap::value_t::const_iterator end = elements.end();
    std::set<IndexMap::value_t::value_type> elementsWithSurfacePhysics;
    std::set<SurfacePhysicsHandler*> surfacePhysicsHandlers;
    std::pair<double, double> currentRange(0.0, 0.0);

    while (elements.size() > 0) {
        auto it = elements.begin();
        if ((*it)->hasSurfacePhysics()) {
            elementsWithSurfacePhysics.insert(*it);
            surfacePhysicsHandlers.insert((*it)->getSurfacePhysics());

            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());
781
        }
782 783

        elements.erase(it);
784
    }
785

786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820
    if (elementsWithSurfacePhysics.size() > 0) {
        std::set<SurfacePhysicsHandler*> oldSPHandlers;
        std::vector<SurfacePhysicsHandler*> leftBehindSPHandlers, newSPHandlers;
        for (auto it: activeSurfacePhysicsHandlers_m) {
            oldSPHandlers.insert(it);
        }

        leftBehindSPHandlers.resize(std::max(oldSPHandlers.size(),
                                             surfacePhysicsHandlers.size()));
        auto last = std::set_difference(oldSPHandlers.begin(), oldSPHandlers.end(),
                                        surfacePhysicsHandlers.begin(), surfacePhysicsHandlers.end(),
                                        leftBehindSPHandlers.begin());
        leftBehindSPHandlers.resize(last - leftBehindSPHandlers.begin());

        for (auto it: leftBehindSPHandlers) {
            if (!it->stillActive()) {
                activeSurfacePhysicsHandlers_m.erase(it);
            }
        }

        newSPHandlers.resize(std::max(oldSPHandlers.size(),
                                      elementsWithSurfacePhysics.size()));
        last = std::set_difference(surfacePhysicsHandlers.begin(), surfacePhysicsHandlers.end(),
                                   oldSPHandlers.begin(), oldSPHandlers.end(),
                                   newSPHandlers.begin());
        newSPHandlers.resize(last - newSPHandlers.begin());

        for (auto it: newSPHandlers) {
            activeSurfacePhysicsHandlers_m.insert(it);
        }

        if(!surfaceStatus_m) {
            msg << level2 << "============== START SURFACE PHYSICS CALCULATION =============" << endl;
            surfaceStatus_m = true;
        }
821
    }
822

823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848
    if (surfaceStatus_m) {
        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
            SurfacePhysicsHandler* onlyDegraderWithParticles = NULL;
            int degradersWithParticlesCount = 0;
            for (auto it: activeSurfacePhysicsHandlers_m) {
                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
            int maxPerNode = itsBunch_m->getLocalNum();
            reduce(maxPerNode, maxPerNode, OpMaxAssign());
            bool allParticlesInMat = (maxPerNode == 0 &&
                                      degradersWithParticlesCount == 1);

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

849
            auto boundingSphere = itsBunch_m->getLocalBoundingSphere();
850 851
            unsigned redifusedParticles = 0;
            for (auto it: activeSurfacePhysicsHandlers_m) {
852 853 854 855 856 857 858 859 860 861
                ElementBase* element = it->getElement();
                CoordinateSystemTrafo refToLocalCSTrafo = (element->getMisalignment() *
                                                           (element->getCSTrafoGlobal2Local() * referenceToLabCSTrafo_m));
                CoordinateSystemTrafo localToRefCSTrafo = refToLocalCSTrafo.inverted();

                const unsigned int localNum = itsBunch_m->getLocalNum();
                for (unsigned int i = 0; i < localNum; ++i) {
                    itsBunch_m->R[i] = refToLocalCSTrafo.transformTo(itsBunch_m->R[i]);
                    itsBunch_m->P[i] = refToLocalCSTrafo.rotateTo(itsBunch_m->P[i]);
                }
862
                boundingSphere.first = refToLocalCSTrafo.transformTo(boundingSphere.first);
863

864
                it->apply(*itsBunch_m, boundingSphere, totalParticlesInSimulation_m);
865 866
                it->print(msg);