ParallelTTracker.h 17.2 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 20 21
#ifndef OPAL_ParallelTTracker_HH
#define OPAL_ParallelTTracker_HH

// ------------------------------------------------------------------------
// $RCSfile: ParallelTTracker.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelTTracker
//
// ------------------------------------------------------------------------
//
// $Date: 2004/11/12 20:10:11 $
// $Author: adelmann $
//
// ------------------------------------------------------------------------

#include "Algorithms/Tracker.h"
22
#include "Steppers/BorisPusher.h"
gsell's avatar
gsell committed
23 24
#include "Structure/DataSink.h"

25
#include "BasicActions/Option.h"
26
#include "Utilities/Options.h"
27

gsell's avatar
gsell committed
28 29
#include "Physics/Physics.h"

30 31
#include "Algorithms/OrbitThreader.h"
#include "Algorithms/IndexMap.h"
32 33
#include "AbsBeamline/AlignWrapper.h"
#include "AbsBeamline/BeamBeam.h"
34
#include "AbsBeamline/Collimator.h"
35 36
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Diagnostic.h"
37
#include "AbsBeamline/Degrader.h"
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/ElementBase.h"
#include "AbsBeamline/Lambertson.h"
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RBend3D.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/TravelingWave.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/ParallelPlate.h"
#include "AbsBeamline/CyclotronValley.h"
gsell's avatar
gsell committed
56 57

#include "Beamlines/Beamline.h"
58
#include "Beamlines/FlaggedBeamline.h"
gsell's avatar
gsell committed
59 60 61 62 63
#include "Elements/OpalBeamline.h"
#include "Solvers/WakeFunction.hh"

#include <list>
#include <vector>
64
#include <queue>
gsell's avatar
gsell committed
65

66
class BorisPusher;
67
class ParticleMatterInteractionHandler;
gsell's avatar
gsell committed
68 69 70 71 72 73 74 75 76 77

class ParallelTTracker: public Tracker {

public:
    /// Constructor.
    //  The beam line to be tracked is "bl".
    //  The particle reference data are taken from "data".
    //  The particle bunch tracked is initially empty.
    //  If [b]revBeam[/b] is true, the beam runs from s = C to s = 0.
    //  If [b]revTrack[/b] is true, we track against the beam.
78 79 80 81
    explicit ParallelTTracker(const Beamline &bl,
                              const PartData &data,
                              bool revBeam,
                              bool revTrack);
gsell's avatar
gsell committed
82 83 84 85 86 87 88

    /// Constructor.
    //  The beam line to be tracked is "bl".
    //  The particle reference data are taken from "data".
    //  The particle bunch tracked is taken from [b]bunch[/b].
    //  If [b]revBeam[/b] is true, the beam runs from s = C to s = 0.
    //  If [b]revTrack[/b] is true, we track against the beam.
89
    explicit ParallelTTracker(const Beamline &bl,
frey_m's avatar
frey_m committed
90
                              PartBunch *bunch,
91 92 93 94 95 96 97 98 99
                              DataSink &ds,
                              const PartData &data,
                              bool revBeam,
                              bool revTrack,
                              const std::vector<unsigned long long> &maxSTEPS,
                              double zstart,
                              const std::vector<double> &zstop,
                              const std::vector<double> &dt);

gsell's avatar
gsell committed
100 101 102 103 104 105 106 107 108 109 110

    virtual ~ParallelTTracker();

    virtual void visitAlignWrapper(const AlignWrapper &);

    /// Apply the algorithm to a BeamBeam.
    virtual void visitBeamBeam(const BeamBeam &);

    /// Apply the algorithm to a collimator.
    virtual void visitCollimator(const Collimator &);

adelmann's avatar
adelmann committed
111

gsell's avatar
gsell committed
112 113 114
    /// Apply the algorithm to a Corrector.
    virtual void visitCorrector(const Corrector &);

adelmann's avatar
adelmann committed
115 116 117
    /// Apply the algorithm to a Degrader.
    virtual void visitDegrader(const Degrader &);

gsell's avatar
gsell committed
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
    /// Apply the algorithm to a Diagnostic.
    virtual void visitDiagnostic(const Diagnostic &);

    /// Apply the algorithm to a Drift.
    virtual void visitDrift(const Drift &);

    /// Apply the algorithm to a Lambertson.
    virtual void visitLambertson(const Lambertson &);

    /// Apply the algorithm to a Marker.
    virtual void visitMarker(const Marker &);

    /// Apply the algorithm to a Monitor.
    virtual void visitMonitor(const Monitor &);

    /// Apply the algorithm to a Multipole.
    virtual void visitMultipole(const Multipole &);

    /// Apply the algorithm to a Probe.
    virtual void visitProbe(const Probe &);

    /// Apply the algorithm to a RBend.
    virtual void visitRBend(const RBend &);

142 143 144
    /// Apply the algorithm to a RBend.
    virtual void visitRBend3D(const RBend3D &);

gsell's avatar
gsell committed
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
    /// Apply the algorithm to a RFCavity.
    virtual void visitRFCavity(const RFCavity &);

    /// Apply the algorithm to a RFCavity.
    virtual void visitTravelingWave(const TravelingWave &);

    /// Apply the algorithm to a RFQuadrupole.
    virtual void visitRFQuadrupole(const RFQuadrupole &);

    /// Apply the algorithm to a SBend.
    virtual void visitSBend(const SBend &);

    /// Apply the algorithm to a Separator.
    virtual void visitSeparator(const Separator &);

    /// Apply the algorithm to a Septum.
    virtual void visitSeptum(const Septum &);

    /// Apply the algorithm to a Solenoid.
    virtual void visitSolenoid(const Solenoid &);

166 167 168
    /// Apply the algorithm to a Solenoid.
    virtual void visitSource(const Source &);

gsell's avatar
gsell committed
169 170 171 172 173 174 175 176 177 178 179
    /// Apply the algorithm to a ParallelPlate.
    virtual void visitParallelPlate(const ParallelPlate &);

    /// Apply the algorithm to a CyclotronValley.
    virtual void visitCyclotronValley(const CyclotronValley &);


    /// Apply the algorithm to a beam line.
    //  overwrite the execute-methode from DefaultVisitor
    virtual void visitBeamline(const Beamline &);

180 181 182
    /// Apply the algorithm to the top-level beamline.
    //  overwrite the execute-methode from DefaultVisitor
    virtual void execute();
gsell's avatar
gsell committed
183 184 185 186 187 188 189 190

private:

    // Not implemented.
    ParallelTTracker();
    ParallelTTracker(const ParallelTTracker &);
    void operator=(const ParallelTTracker &);

191 192
    /******************** STATE VARIABLES ***********************************/

193
    DataSink *itsDataSink_m;
gsell's avatar
gsell committed
194 195

    OpalBeamline itsOpalBeamline_m;
196

197 198
    Vector_t RefPartR_m;
    Vector_t RefPartP_m;
199

200 201 202 203
    bool globalEOL_m;

    bool wakeStatus_m;

204
    bool deletedParticles_m;
205

206
    WakeFunction* wakeFunction_m;
207

208
    double pathLength_m;
209

210 211
    /// where to start
    double zstart_m;
212 213

    /// where to stop
214
    std::queue<double> zStop_m;
215

216
    double dtCurrentTrack_m;
217
    std::queue<double> dtAllTracks_m;
218

219 220
    /// The maximal number of steps the system is integrated per TRACK
    std::queue<unsigned long long> localTrackSteps_m;
adelmann's avatar
adelmann committed
221

222 223 224 225 226 227 228 229 230 231
    // This variable controls the minimal number of steps of emission (using bins)
    // before we can merge the bins
    int minStepforReBin_m;

    // The space charge solver crashes if we use less than ~10 particles.
    // This variable controls the number of particles to be emitted before we use
    // the space charge solver.
    size_t minBinEmitted_m;

    // this variable controls the minimal number of steps until we repartition the particles
232
    unsigned int repartFreq_m;
233 234 235 236 237 238 239

    unsigned int emissionSteps_m;

    size_t numParticlesInSimulation_m;

    IpplTimings::TimerRef timeIntegrationTimer1_m;
    IpplTimings::TimerRef timeIntegrationTimer2_m;
240
    IpplTimings::TimerRef fieldEvaluationTimer_m ;
241 242
    IpplTimings::TimerRef BinRepartTimer_m;
    IpplTimings::TimerRef WakeFieldTimer_m;
243

244
    CoordinateSystemTrafo referenceToLabCSTrafo_m;
245

246
    std::set<ParticleMatterInteractionHandler*> activeParticleMatterInteractionHandlers_m;
247
    bool particleMaterStatus_m;
248 249 250

    unsigned long totalParticlesInSimulation_m;

251 252 253 254 255 256 257 258 259 260 261 262
#ifdef OPAL_DKS
    DKSOPAL *dksbase;
    void *r_ptr;
    void *p_ptr;
    void *Ef_ptr;
    void *Bf_ptr;
    void *dt_ptr;

    int stream1;
    int stream2;

    int numDeviceElements;
263

264 265 266 267 268 269 270 271 272 273
    void setupDKS();
    void allocateDeviceMemory();
    void freeDeviceMemory();
    void sendRPdt();
    void sendEfBf();
    void pushParticlesDKS(bool send = true);
    void kickParticlesDKS();

#endif

274 275
    /********************** END VARIABLES ***********************************/

gsell's avatar
gsell committed
276
    void kickParticles(const BorisPusher &pusher);
277 278
    void pushParticles(const BorisPusher &pusher);
    void updateReferenceParticle(const BorisPusher &pusher);
gsell's avatar
gsell committed
279

280
    void writePhaseSpace(const long long step, bool psDump, bool statDump);
gsell's avatar
gsell committed
281

282
    /********** BEGIN AUTOPHSING STUFF **********/
283
    void updateRFElement(std::string elName, double maxPhi);
284
    void printRFPhases();
285 286
    void saveCavityPhases();
    void restoreCavityPhases();
287 288
    /************ END AUTOPHSING STUFF **********/

289
    void prepareSections();
290

291 292
    void timeIntegration1(BorisPusher & pusher);
    void timeIntegration2(BorisPusher & pusher);
293
    void selectDT();
294
    void changeDT();
295
    void emitParticles(long long step);
296
    void computeExternalFields(OrbitThreader &oth);
297 298
    void computeWakefield(IndexMap::value_t &elements);
    void computeParticleMatterInteraction(IndexMap::value_t elements, OrbitThreader &oth);
299
    void computeSpaceChargeFields(unsigned long long step);
300
    // void prepareOpalBeamlineSections();
301
    void dumpStats(long long step, bool psDump, bool statDump);
302 303 304 305 306
    void setOptionalVariables();
    bool hasEndOfLineReached();
    void handleRestartRun();
    void prepareEmission();
    void setTime();
307
    void doBinaryRepartition();
gsell's avatar
gsell committed
308

309
    void transformBunch(const CoordinateSystemTrafo &trafo);
310

311 312 313
    void updateRefToLabCSTrafo(const BorisPusher &pusher);
    void findStartPosition(const BorisPusher &pusher);
    void autophaseCavities(const BorisPusher &pusher);
gsell's avatar
gsell committed
314

315
    void evenlyDistributeParticles();
gsell's avatar
gsell committed
316

317 318 319 320 321
    static unsigned long long getMaxSteps(std::queue<unsigned long long> numSteps);

    // std::ofstream logger_m;
    // size_t loggingFrequency_m;
};
gsell's avatar
gsell committed
322 323

inline void ParallelTTracker::visitAlignWrapper(const AlignWrapper &wrap) {
324
    itsOpalBeamline_m.visit(wrap, *this, itsBunch_m);
gsell's avatar
gsell committed
325 326 327
}

inline void ParallelTTracker::visitBeamBeam(const BeamBeam &bb) {
328
    itsOpalBeamline_m.visit(bb, *this, itsBunch_m);
gsell's avatar
gsell committed
329 330 331 332
}


inline void ParallelTTracker::visitCollimator(const Collimator &coll) {
333
    itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
gsell's avatar
gsell committed
334 335 336 337
}


inline void ParallelTTracker::visitCorrector(const Corrector &corr) {
338
    itsOpalBeamline_m.visit(corr, *this, itsBunch_m);
gsell's avatar
gsell committed
339 340 341
}


adelmann's avatar
adelmann committed
342
inline void ParallelTTracker::visitDegrader(const Degrader &deg) {
343
    itsOpalBeamline_m.visit(deg, *this, itsBunch_m);
adelmann's avatar
adelmann committed
344 345 346
}


gsell's avatar
gsell committed
347
inline void ParallelTTracker::visitDiagnostic(const Diagnostic &diag) {
348
    itsOpalBeamline_m.visit(diag, *this, itsBunch_m);
gsell's avatar
gsell committed
349 350 351 352
}


inline void ParallelTTracker::visitDrift(const Drift &drift) {
353
    itsOpalBeamline_m.visit(drift, *this, itsBunch_m);
gsell's avatar
gsell committed
354 355 356 357
}


inline void ParallelTTracker::visitLambertson(const Lambertson &lamb) {
358
    itsOpalBeamline_m.visit(lamb, *this, itsBunch_m);
gsell's avatar
gsell committed
359 360 361 362
}


inline void ParallelTTracker::visitMarker(const Marker &marker) {
363
    itsOpalBeamline_m.visit(marker, *this, itsBunch_m);
gsell's avatar
gsell committed
364 365 366 367
}


inline void ParallelTTracker::visitMonitor(const Monitor &mon) {
368
    itsOpalBeamline_m.visit(mon, *this, itsBunch_m);
gsell's avatar
gsell committed
369 370 371 372
}


inline void ParallelTTracker::visitMultipole(const Multipole &mult) {
373
    itsOpalBeamline_m.visit(mult, *this, itsBunch_m);
gsell's avatar
gsell committed
374 375 376
}

inline void ParallelTTracker::visitProbe(const Probe &prob) {
377
    itsOpalBeamline_m.visit(prob, *this, itsBunch_m);
gsell's avatar
gsell committed
378 379 380 381
}


inline void ParallelTTracker::visitRBend(const RBend &bend) {
382 383 384 385 386
    itsOpalBeamline_m.visit(bend, *this, itsBunch_m);
}

inline void ParallelTTracker::visitRBend3D(const RBend3D &bend) {
    itsOpalBeamline_m.visit(bend, *this, itsBunch_m);
gsell's avatar
gsell committed
387 388 389 390
}


inline void ParallelTTracker::visitRFCavity(const RFCavity &as) {
391
    itsOpalBeamline_m.visit(as, *this, itsBunch_m);
gsell's avatar
gsell committed
392 393 394
}

inline void ParallelTTracker::visitTravelingWave(const TravelingWave &as) {
395
    itsOpalBeamline_m.visit(as, *this, itsBunch_m);
gsell's avatar
gsell committed
396 397 398 399
}


inline void ParallelTTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
400
    itsOpalBeamline_m.visit(rfq, *this, itsBunch_m);
gsell's avatar
gsell committed
401 402 403
}

inline void ParallelTTracker::visitSBend(const SBend &bend) {
404
    itsOpalBeamline_m.visit(bend, *this, itsBunch_m);
gsell's avatar
gsell committed
405 406 407 408
}


inline void ParallelTTracker::visitSeparator(const Separator &sep) {
409
    itsOpalBeamline_m.visit(sep, *this, itsBunch_m);
gsell's avatar
gsell committed
410 411 412 413
}


inline void ParallelTTracker::visitSeptum(const Septum &sept) {
414
    itsOpalBeamline_m.visit(sept, *this, itsBunch_m);
gsell's avatar
gsell committed
415 416 417 418
}


inline void ParallelTTracker::visitSolenoid(const Solenoid &solenoid) {
419
    itsOpalBeamline_m.visit(solenoid, *this, itsBunch_m);
gsell's avatar
gsell committed
420 421
}

422 423 424
inline void ParallelTTracker::visitSource(const Source &source) {
    itsOpalBeamline_m.visit(source, *this, itsBunch_m);
}
gsell's avatar
gsell committed
425 426

inline void ParallelTTracker::visitParallelPlate(const ParallelPlate &pplate) {
427
    itsOpalBeamline_m.visit(pplate, *this, itsBunch_m);
gsell's avatar
gsell committed
428 429 430
}

inline void ParallelTTracker::visitCyclotronValley(const CyclotronValley &cv) {
431
    itsOpalBeamline_m.visit(cv, *this, itsBunch_m);
gsell's avatar
gsell committed
432 433 434
}

inline void ParallelTTracker::kickParticles(const BorisPusher &pusher) {
435 436 437
    int localNum = itsBunch_m->getLocalNum();
    for (int i = 0; i < localNum; ++i)
        pusher.kick(itsBunch_m->R[i], itsBunch_m->P[i], itsBunch_m->Ef[i], itsBunch_m->Bf[i], itsBunch_m->dt[i]);
gsell's avatar
gsell committed
438 439
}

440 441
inline void ParallelTTracker::pushParticles(const BorisPusher &pusher) {
    itsBunch_m->switchToUnitlessPositions(true);
gsell's avatar
gsell committed
442

443 444
    for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
        pusher.push(itsBunch_m->R[i], itsBunch_m->P[i], itsBunch_m->dt[i]);
gsell's avatar
gsell committed
445
    }
446
    itsBunch_m->switchOffUnitlessPositions(true);
gsell's avatar
gsell committed
447 448
}

449 450 451
inline
unsigned long long ParallelTTracker::getMaxSteps(std::queue<unsigned long long> numSteps) {
    unsigned long long totalNumSteps = 0;
gsell's avatar
gsell committed
452

453 454 455 456
    while (numSteps.size() > 0) {
        totalNumSteps += numSteps.front();
        numSteps.pop();
    }
gsell's avatar
gsell committed
457

458
    return totalNumSteps;
gsell's avatar
gsell committed
459 460
}

461 462
#ifdef OPAL_DKS
inline void ParallelTTracker::setupDKS() {
463 464 465 466
    dksbase = new DKSOPAL;
    dksbase->setAPI("Cuda");
    dksbase->setDevice("-gpu");
    dksbase->initDevice();
467

468 469
    dksbase->createStream(stream1);
    dksbase->createStream(stream2);
470 471 472 473
}


inline void ParallelTTracker::allocateDeviceMemory() {
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491

    int ierr;
    r_ptr = dksbase->allocateMemory<Vector_t>(itsBunch_m->getLocalNum(), ierr);
    p_ptr = dksbase->allocateMemory<Vector_t>(itsBunch_m->getLocalNum(), ierr);
    Ef_ptr = dksbase->allocateMemory<Vector_t>(itsBunch_m->getLocalNum(), ierr);
    Bf_ptr = dksbase->allocateMemory<Vector_t>(itsBunch_m->getLocalNum(), ierr);
    dt_ptr = dksbase->allocateMemory<double>(itsBunch_m->getLocalNum(), ierr);

    if (Ippl::getNodes() == 1) {
        dksbase->registerHostMemory(&itsBunch_m->R[0], itsBunch_m->getLocalNum());
        dksbase->registerHostMemory(&itsBunch_m->P[0], itsBunch_m->getLocalNum());
        dksbase->registerHostMemory(&itsBunch_m->Ef[0], itsBunch_m->getLocalNum());
        dksbase->registerHostMemory(&itsBunch_m->Bf[0], itsBunch_m->getLocalNum());
        dksbase->registerHostMemory(&itsBunch_m->dt[0], itsBunch_m->getLocalNum());

    }

    numDeviceElements = itsBunch_m->getLocalNum();
492 493 494
}

inline  void ParallelTTracker::freeDeviceMemory() {
495 496 497 498 499 500 501 502 503 504 505 506 507
    dksbase->freeMemory<Vector_t>(r_ptr, numDeviceElements);
    dksbase->freeMemory<Vector_t>(p_ptr, numDeviceElements);
    dksbase->freeMemory<Vector_t>(Ef_ptr, numDeviceElements);
    dksbase->freeMemory<Vector_t>(Bf_ptr, numDeviceElements);
    dksbase->freeMemory<double>(dt_ptr, numDeviceElements);

    if (Ippl::getNodes() == 1) {
        dksbase->unregisterHostMemory(&itsBunch_m->R[0]);
        dksbase->unregisterHostMemory(&itsBunch_m->P[0]);
        dksbase->unregisterHostMemory(&itsBunch_m->Ef[0]);
        dksbase->unregisterHostMemory(&itsBunch_m->Bf[0]);
        dksbase->unregisterHostMemory(&itsBunch_m->dt[0]);
    }
508 509 510
}

inline void ParallelTTracker::sendRPdt() {
511 512 513
    dksbase->writeDataAsync<Vector_t>(r_ptr, &itsBunch_m->R[0], itsBunch_m->getLocalNum(), stream1);
    dksbase->writeDataAsync<Vector_t>(p_ptr, &itsBunch_m->P[0], itsBunch_m->getLocalNum(), stream1);
    dksbase->writeDataAsync<double>(dt_ptr, &itsBunch_m->dt[0], itsBunch_m->getLocalNum(), stream1);
514 515 516
}

inline void ParallelTTracker::sendEfBf() {
517 518 519 520
    dksbase->writeDataAsync<Vector_t>(Ef_ptr, &itsBunch_m->Ef[0],
                                      itsBunch_m->getLocalNum(), stream1);
    dksbase->writeDataAsync<Vector_t>(Bf_ptr, &itsBunch_m->Bf[0],
                                      itsBunch_m->getLocalNum(), stream1);
521 522 523 524
}

inline void ParallelTTracker::pushParticlesDKS(bool send) {

525 526 527 528 529 530 531 532
    //send data to the GPU
    if (send)
        sendRPdt();
    //execute particle push on the GPU
    dksbase->callParallelTTrackerPush(r_ptr, p_ptr, dt_ptr, itsBunch_m->getLocalNum(),
                                      Physics::c, stream1);
    //get particles back to CPU
    dksbase->readData<Vector_t>(r_ptr, &itsBunch_m->R[0], itsBunch_m->getLocalNum(), stream1);
533 534 535 536
}

inline
void ParallelTTracker::kickParticlesDKS() {
537 538 539 540 541 542
    //send data to the GPU
    sendEfBf();
    sendRPdt();

    //sync device
    dksbase->syncDevice();
543

544 545 546 547
    //execute kick on the GPU
    dksbase->callParallelTTrackerKick(r_ptr, p_ptr, Ef_ptr, Bf_ptr, dt_ptr,
                                      itsReference.getQ(), itsReference.getM(),
                                      itsBunch_m->getLocalNum(), Physics::c, stream2);
548

549
    dksbase->syncDevice();
550

551 552
    //get data back from GPU
    dksbase->readDataAsync<Vector_t>(p_ptr, &itsBunch_m->P[0], itsBunch_m->getLocalNum(), stream2);
553 554 555 556

}
#endif

557
#endif // OPAL_ParallelTTracker_HH