ParallelTracker.h

Go to the documentation of this file.

//
// Class ParallelTracker
//   OPAL-T tracker.
//   The visitor class for tracking particles with time as independent
//   variable.
//
// Copyright (c) 200x - 2014, Christof Kraus, Paul Scherrer Institut, Villigen PSI, Switzerland
//               2015 - 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
//               2017 - 2020, Christof Metzger-Kraus
// All rights reserved
//
// 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/>.
//
#ifndef OPAL_ParallelTracker_HH
#define OPAL_ParallelTracker_HH
 
#include "Algorithms/StepSizeConfig.h"
#include "Algorithms/Tracker.h"
#include "Steppers/BorisPusher.h"
#include "Structure/DataSink.h"
 
#include "BasicActions/Option.h"
#include "Utilities/Options.h"
 
#include "Physics/Physics.h"
 
#include "Algorithms/IndexMap.h"
#include "Algorithms/OrbitThreader.h"
 
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/ElementBase.h"
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Multipole.h"
#include "AbsBeamline/MultipoleT.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/Ring.h"
#include "AbsBeamline/ScalingFFAMagnet.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/TravelingWave.h"
#include "Beamlines/Beamline.h"
#include "Elements/OpalBeamline.h"
 
#include <list>
#include <memory>
#include <tuple>
#include <vector>
 
class ParticleMatterInteractionHandler;
class PluginElement;
 
class ParallelTracker : public Tracker {
 
    DataSink* itsDataSink_m;
 
    OpalBeamline itsOpalBeamline_m;
 
    /*
      Ring specifics
    */
 
    // we store a pointer explicitly to the Ring
    Ring* opalRing_m;
 
    bool globalEOL_m;
 
    bool wakeStatus_m;
 
    bool deletedParticles_m;
 
    WakeFunction* wakeFunction_m;
 
    double pathLength_m;

    double zstart_m;

    StepSizeConfig stepSizes_m;
 
    double dtCurrentTrack_m;
 
    // This variable controls the minimal number of steps of emission (using bins)
    // before we can merge the bins
    int minStepforReBin_m;
 
    // this variable controls the minimal number of steps until we repartition the particles
    unsigned long long repartFreq_m;
 
    unsigned int emissionSteps_m;
 
    size_t numParticlesInSimulation_m;
 
    IpplTimings::TimerRef timeIntegrationTimer1_m;
    IpplTimings::TimerRef timeIntegrationTimer2_m;
    IpplTimings::TimerRef fieldEvaluationTimer_m;
    IpplTimings::TimerRef WakeFieldTimer_m;
    IpplTimings::TimerRef PluginElemTimer_m;
    IpplTimings::TimerRef BinRepartTimer_m;
    IpplTimings::TimerRef OrbThreader_m;
    
    std::set<ParticleMatterInteractionHandler*> activeParticleMatterInteractionHandlers_m;
    bool particleMatterStatus_m;
 
public:
    typedef std::vector<double> dvector_t;
    typedef std::vector<int> ivector_t;
    typedef std::pair<double[8], Component*> element_pair;
    typedef std::pair<ElementType, element_pair> type_pair;
    typedef std::list<type_pair*> beamline_list;

    //  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.
    explicit ParallelTracker(const Beamline& bl, const PartData& data, bool revBeam, bool revTrack);

    //  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.
    explicit ParallelTracker(
        const Beamline& bl, PartBunch_t* bunch, 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);
 
    virtual ~ParallelTracker();

    //  overwrite the execute-methode from DefaultVisitor
    virtual void execute();

    //  overwrite the execute-methode from DefaultVisitor
    virtual void visitBeamline(const Beamline&);

    virtual void visitDrift(const Drift&);

    virtual void visitRing(const Ring& ring);

    virtual void visitMarker(const Marker&);

    virtual void visitMultipole(const Multipole&);

    virtual void visitMultipoleT(const MultipoleT&);

    virtual void visitOffset(const Offset&);

    virtual void visitRFCavity(const RFCavity&);

    virtual void visitSolenoid(const Solenoid&);

    virtual void visitTravelingWave(const TravelingWave&);

    virtual void visitScalingFFAMagnet(const ScalingFFAMagnet& bend);

    virtual void visitVerticalFFAMagnet(const VerticalFFAMagnet& bend);
 
    // made following public: __host__ __device__ lambda cannot have private or protected access within its class
    void kickParticles(const BorisPusher& pusher);
 
    void pushParticles(const BorisPusher& pusher);
 
    void timeIntegration2(BorisPusher& pusher);
 
    void changeDT(bool backTrack = false);
 
    void computeSpaceChargeFields(unsigned long long step);
 
    void setTime();
private:
    // Not implemented.
    ParallelTracker();
    ParallelTracker(const ParallelTracker&);
    void operator=(const ParallelTracker&);
 
    /*
      Ring specifics
    */
 
    unsigned turnnumber_m;
 
    std::list<Component*> myElements;
    beamline_list FieldDimensions;
    double bega;
    double referenceR;
    double referenceTheta;
    double referenceZ = 0.0;
 
    double referencePr;
    double referencePt;
    double referencePz = 0.0;
    double referencePtot;
 
    double referencePsi;
    double referencePhi;
 
    double sinRefTheta_m;
    double cosRefTheta_m;
 
    void buildupFieldList(double BcParameter[], ElementType elementType, Component* elptr);
 
    std::vector<PluginElement*> pluginElements_m;
 
 
 
 
 
    /********************** END VARIABLES ***********************************/
 
    void updateReferenceParticle(const BorisPusher& pusher);
 
    void writePhaseSpace(const long long step, bool psDump, bool statDump);
 
    /********** BEGIN AUTOPHSING STUFF **********/
    void updateRFElement(std::string elName, double maxPhi);
    void printRFPhases();
    void saveCavityPhases();
    void restoreCavityPhases();
    /************ END AUTOPHSING STUFF **********/
 
    void prepareSections();
 
    void timeIntegration1(BorisPusher& pusher);
    void selectDT(bool backTrack = false);
    void emitParticles(long long step);
    void computeExternalFields(OrbitThreader& oth);
    void computeWakefield(IndexMap::value_t& elements);
    void computeParticleMatterInteraction(IndexMap::value_t elements, OrbitThreader& oth);
 
    // void prepareOpalBeamlineSections();
    void dumpStats(long long step, bool psDump, bool statDump);
    void setOptionalVariables();
    bool hasEndOfLineReached(const BoundingBox& globalBoundingBox);
    void handleRestartRun();
 
    void doBinaryRepartition();
 
    void transformBunch(const CoordinateSystemTrafo& trafo);
 
    void updateReference(const BorisPusher& pusher);
    void updateRefToLabCSTrafo();
    void applyFractionalStep(const BorisPusher& pusher, double tau);
    void findStartPosition(const BorisPusher& pusher);
    void autophaseCavities(const BorisPusher& pusher);
 
    bool applyPluginElements(const double dt);
};
 
inline void ParallelTracker::visitDrift(const Drift& drift) {
    itsOpalBeamline_m.visit(drift, *this, itsBunch_m);
}
 
inline void ParallelTracker::visitMarker(const Marker& marker) {
    itsOpalBeamline_m.visit(marker, *this, itsBunch_m);
}
 
inline void ParallelTracker::visitMultipole(const Multipole& mult) {
    itsOpalBeamline_m.visit(mult, *this, itsBunch_m);
}
 
inline void ParallelTracker::visitMultipoleT(const MultipoleT& mult) {
    itsOpalBeamline_m.visit(mult, *this, itsBunch_m);
}
 
inline void ParallelTracker::visitRFCavity(const RFCavity& as) {
    itsOpalBeamline_m.visit(as, *this, itsBunch_m);
}
 
inline void ParallelTracker::visitSolenoid(const Solenoid& so) {
    itsOpalBeamline_m.visit(so, *this, itsBunch_m);
}
 
inline void ParallelTracker::visitTravelingWave(const TravelingWave& as) {
    itsOpalBeamline_m.visit(as, *this, itsBunch_m);
}
 
inline void ParallelTracker::kickParticles(const BorisPusher& pusher) {
    
    auto Rview  = itsBunch_m->getParticleContainer()->R.getView();
    auto Pview  = itsBunch_m->getParticleContainer()->P.getView();
    auto dtview = itsBunch_m->getParticleContainer()->dt.getView();
    auto Efview = itsBunch_m->getParticleContainer()->E.getView();
    auto Bfview = itsBunch_m->getParticleContainer()->B.getView();
 
 
    const double mass = itsReference.getM();
    const double charge = itsReference.getQ();
 
    Kokkos::parallel_for(
                         "kickParticles", ippl::getRangePolicy(Rview),
                         KOKKOS_LAMBDA(const int i) {
                             Vector_t<double, 3> p = {Pview(i)[0],Pview(i)[1],Pview(i)[2]};
 
                             Vector_t<double, 3> e = {Efview(i)[0],Efview(i)[1],Efview(i)[2]};
                             Vector_t<double, 3> b = {Bfview(i)[0],Bfview(i)[1],Bfview(i)[2]};
                             double dt = dtview(i);
 
                             // pusher.kick(x,p,e,b,dt);
                             // Implementation follows chapter 4-4, p. 61 - 63 from
                             // Birdsall, C. K. and Langdon, A. B. (1985). Plasma physics
                             // via computer simulation.
                             //
                             // Up to finite precision effects, the new implementation is equivalent to the
                             // old one, but uses less floating point operations.
                             //
                             // Relativistic variant implemented below is described in
                             // chapter 15-4, p. 356 - 357.
                             // However, since other units are used here, small
                             // modifications are required. The relativistic variant can be derived
                             // from the nonrelativistic one by replacing
                             //     mass
                             // by
                             //     gamma * rest mass
                             // and transforming the units.
                             //
                             // Parameters:
                             //     R = x / (c * dt): Scaled position x, not used in here
                             //     P = v / c * gamma: Scaled velocity v
                             //     Ef: Electric field
                             //     Bf: Magnetic field
                             //     dt: Timestep
                             //     mass = rest energy = rest mass * c * c
                             //     charge
                             
                             // Half step E
                             p += 0.5 * dt * charge * Physics::c / mass * e;
 
                             // Full step B
 
                             const double gamma          = Kokkos::sqrt(1.0 + dot(p, p));
                             Vector_t<double, 3> const t = 0.5 * dt * charge * Physics::c * Physics::c / (gamma * mass) * b;
                             Vector_t<double, 3> const w = p + cross(p, t);
                             Vector_t<double, 3> const s = 2.0 / (1.0 + dot(t, t)) * t;
                             p += cross(w, s);
 
 
                             // Half step E
                             p += 0.5 * dt * charge * Physics::c / mass * e;
 
                             Pview(i) = p;
                                 
                         });
 
    itsBunch_m->getParticleContainer()->update();
    ippl::Comm->barrier();
}
 
inline void ParallelTracker::pushParticles(const BorisPusher& pusher) {
 
    itsBunch_m->switchToUnitlessPositions(true);
 
    auto Rview  = itsBunch_m->getParticleContainer()->R.getView();
    auto Pview  = itsBunch_m->getParticleContainer()->P.getView();
    auto dtview = itsBunch_m->getParticleContainer()->dt.getView();
 
    Kokkos::parallel_for(
                         "pushParticles", ippl::getRangePolicy(Rview),
                         KOKKOS_LAMBDA(const int i) {
                             Vector_t<double, 3> x = {Rview(i)[0],Rview(i)[1],Rview(i)[2]};
                             Vector_t<double, 3> p = {Pview(i)[0],Pview(i)[1],Pview(i)[2]};
                             double dt = dtview(i);

                             // \TODO check +-
                             
                             // pusher.push(x,p,dt);
                             x = 0.5 * dt * p / Kokkos::sqrt(1.0 + dot(p));
                             Rview(i) += x;
                         });
 
 
    itsBunch_m->switchOffUnitlessPositions(true);
    itsBunch_m->getParticleContainer()->update();
    ippl::Comm->barrier();
}
 
#endif  // OPAL_ParallelTracker_HH