ParallelCyclotronTracker: Introduce Stepper class

src/Algorithms/AutophaseTracker.h

@@ -6,7 +6,7 @@
 //
 
 #include "Algorithms/PartBunch.h"
-#include "Algorithms/PartPusher.h"
+#include "Steppers/BorisPusher.h"
 #include "Algorithms/PartData.h"
 #include "Algorithms/PBunchDefs.h"
 

src/Algorithms/OrbitThreader.h

@@ -4,7 +4,7 @@
 #include "Algorithms/IndexMap.h"
 #include "Algorithms/Vektor.h"
 #include "Elements/OpalBeamline.h"
-#include "Algorithms/PartPusher.h"
+#include "Steppers/BorisPusher.h"
 
 #include <string>
 #include <map>

src/Algorithms/ParallelCyclotronTracker.h

@@ -23,6 +23,8 @@
 #include "AbsBeamline/ElementBase.h"
 #include <vector>
 
+#include "Steppers/Steppers.h"
+
 class BMultipoleField;
 class PartBunch;
 class PlanarArcGeometry;
@@ -45,6 +47,13 @@ struct CavityCrossData {
 class ParallelCyclotronTracker: public Tracker {
 
 public:
+    
+    enum MODE {
+        UNDEFINED = -1,
+        SINGLE = 0,
+        SEO = 1,
+        BUNCH = 2
+    };
 
     typedef std::pair<double[8], Component *>      element_pair;
     typedef std::pair<ElementBase::ElementType, element_pair>        type_pair;
@@ -251,13 +260,6 @@ private:
 
     double PathLength_m;
 
-    // the name of time integrator
-    // The ID of time integrator
-    // 0 --- RK-4(default)
-    // 1 --- LF-2
-    // 2 --- MTS
-    int  timeIntegrator_m;
-
     void Tracker_LF();
     void Tracker_RK4();
     void Tracker_MTS();
@@ -466,17 +468,43 @@ private:
 
     // we store a pointer explicitly to the Ring
     Ring* opalRing_m;
+    
 
     // If Ring is defined take the harmonic number from Ring; else use
     // cyclotron
     double getHarmonicNumber() const;
     
+    typedef std::function<bool(const double&,
+                               const size_t&,
+                               Vector_t&,
+                               Vector_t&)> function_t;
+    
+    std::unique_ptr< Stepper<function_t> > itsStepper_mp;
+    
+    struct settings {
+        int scSolveFreq;
+        int stepsPerTurn;
+        double azimuth_angle0;
+        double azimuth_angle1;
+        double azimuth_angle2;
+        double deltaTheta;
+        int stepsNextCheck;
+    } setup_m;
+    
+    MODE mode_m;
+    
+    stepper::INTEGRATOR stepper_m;
+    
+    void seoMode_m(double& t, const double dt, bool& dumpEachTurn,
+                   std::vector<double>& Ttime, std::vector<double>& Tdeltr,
+                   std::vector<double>& Tdeltz, std::vector<int>& TturnNumber);
     
-    void seo_m();
+    void singleMode_m(double& t, const double dt, bool& dumpEachTurn, double& oldReferenceTheta);
     
-    void single_mode_m();
+    void bunchMode_m(double& t, const double dt, bool& dumpEachTurn);
     
-    void bunch_mode_m();
+    void gapCrossKick_m(size_t i, double t, double dt,
+                        const Vector_t& Rold, const Vector_t& Pold);
 
 };
 

src/Algorithms/ParallelTTracker.cpp

@@ -28,7 +28,6 @@
 #include <limits>
 #include <cmath>
 
-#include "Algorithms/PartPusher.h"
 #include "Algorithms/OrbitThreader.h"
 #include "Algorithms/CavityAutophaser.h"
 #include "Beamlines/Beamline.h"

src/Algorithms/ParallelTTracker.h

@@ -19,7 +19,7 @@
 // ------------------------------------------------------------------------
 
 #include "Algorithms/Tracker.h"
-#include "Algorithms/PartPusher.h"
+#include "Steppers/BorisPusher.h"
 #include "Structure/DataSink.h"
 
 #include "BasicActions/Option.h"

src/CMakeLists.txt

@@ -79,6 +79,7 @@ add_subdirectory (OpalParser)
 add_subdirectory (Match)
 add_subdirectory (PhysicsActions)
 add_subdirectory (Solvers)
+add_subdirectory (Steppers)
 add_subdirectory (Structure)
 add_subdirectory (Tables)
 add_subdirectory (Track)

src/Classic/AbsBeamline/Bend.cpp

@@ -19,7 +19,6 @@
 // ------------------------------------------------------------------------
 
 #include "AbsBeamline/Bend.h"
-#include "Algorithms/PartPusher.h"
 #include "Algorithms/PartBunch.h"
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "Utilities/Options.h"

src/Classic/AbsBeamline/Bend.h

@@ -22,7 +22,7 @@
 // ------------------------------------------------------------------------
 
 #include "AbsBeamline/BendBase.h"
-#include "Algorithms/PartPusher.h"
+#include "Steppers/BorisPusher.h"
 #include "Utilities/GeneralClassicException.h"
 
 #include "gsl/gsl_spline.h"

src/Classic/AbsBeamline/RBend3D.cpp

@@ -15,7 +15,7 @@
 
 #include "AbsBeamline/RBend3D.h"
 #include "Algorithms/PartBunch.h"
-#include "Algorithms/PartPusher.h"
+#include "Steppers/BorisPusher.h"
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "Fields/Fieldmap.h"
 #include "AbstractObjects/OpalData.h"

src/Classic/AbsBeamline/RFCavity.cpp

@@ -21,7 +21,7 @@
 #include "AbsBeamline/RFCavity.h"
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "Algorithms/PartBunch.h"
-#include "Algorithms/PartPusher.h"
+#include "Steppers/BorisPusher.h"
 #include "Fields/Fieldmap.h"
 #include "Utilities/GeneralClassicException.h"
 #include "Utilities/Util.h"

src/Classic/AbsBeamline/SBend.cpp

@@ -19,7 +19,6 @@
 // ------------------------------------------------------------------------
 
 #include "AbsBeamline/SBend.h"
-#include "Algorithms/PartPusher.h"
 #include "Algorithms/PartBunch.h"
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "Utilities/Options.h"

src/Classic/Algorithms/CMakeLists.txt

@@ -48,7 +48,6 @@ set (HDRS
     PartBunch.h
     PartData.h
     Particle.h
-    PartPusher.h
     PBunchDefs.h
     PolynomialTimeDependence.h
     Quaternion.h

src/Classic/Algorithms/PartPusher.h

-#ifndef CLASSIC_PartPusher_H
-#define CLASSIC_PartPusher_H
-
-#include "Algorithms/Vektor.h"
-#include "Algorithms/PartData.h"
-#include "Physics/Physics.h"
-
-/*
-
-
- */
-
-class BorisPusher {
-public:
-    BorisPusher(const PartData &ref);
-    BorisPusher();
-    void initialise(const PartData *ref);
-    void kick(const Vector_t &R, Vector_t &P, const Vector_t &Ef, const Vector_t &Bf, const double &dt) const;
-    void kick(const Vector_t &R, Vector_t &P, const Vector_t &Ef, const Vector_t &Bf, const double &dt, const double &mass, const double &charge) const;
-    void push(Vector_t &R, const Vector_t &P, const double &dt) const;
-private:
-    const PartData *itsReference;
-};
-
-inline BorisPusher::BorisPusher(const PartData &ref):
-    itsReference(&ref)
-{ }
-
-inline BorisPusher::BorisPusher():
-    itsReference(NULL)
-{}
-
-inline void BorisPusher::initialise(const PartData *ref)
-{ itsReference = ref; }
-
-inline void BorisPusher::kick(const Vector_t &R, Vector_t &P, const Vector_t &Ef, const Vector_t &Bf, const double &dt) const {
-
-	kick(R, P, Ef, Bf, dt, itsReference->getM(), itsReference->getQ());
-
-}
-
-
-inline void BorisPusher::kick(const Vector_t &R, Vector_t &P, const Vector_t &Ef, const Vector_t &Bf, const double &dt, const double &mass, const double &charge) const {
-
-	// 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
-
-	using Physics::c;
-
-	// Half step E
-	P += 0.5 * dt * charge * c / mass * Ef;
-
-	// Full step B
-
-	/*
-             LF
-	double const gamma = sqrt(1.0 + dot(P, P));
-	Vector_t const t = dt * charge * c * c / (gamma * mass) * Bf;
-	P +=  cross(P, t);
-	*/
-
-	double const gamma = sqrt(1.0 + dot(P, P));
-	Vector_t const t = 0.5 * dt * charge * c * c / (gamma * mass) * Bf;
-	Vector_t const w = P + cross(P, t);
-	Vector_t const s = 2.0 / (1.0 + dot(t, t)) * t;
-	P += cross(w, s);
-
-	/* a poor Leap-Frog
-	   P += 1.0 * dt * charge * c * c / (gamma * mass) * cross(P,Bf);
-	*/
-
-	// Half step E
-	P += 0.5 * dt * charge * c / mass * Ef;
-
-}
-
-
-inline void BorisPusher::push(Vector_t &R, const Vector_t &P, const double &/* dt */) const {
-    /** \f[ \vec{x}_{n+1/2} = \vec{x}_{n} + \frac{1}{2}\vec{v}_{n-1/2}\quad (= \vec{x}_{n} + \frac{\Delta t}{2} \frac{\vec{\beta}_{n-1/2}\gamma_{n-1/2}}{\gamma_{n-1/2}}) \f]
-     *
-     * \code
-     * R[i] += 0.5 * P[i] * recpgamma;
-     * \endcode
-     */
-    R += 0.5 * P / sqrt(1.0 + dot(P, P));
-}
-
-#endif
\ No newline at end of file

src/Steppers/BorisPusher.h

+#ifndef CLASSIC_PartPusher_H
+#define CLASSIC_PartPusher_H
+
+#include "Algorithms/Vektor.h"
+#include "Algorithms/PartData.h"
+#include "Physics/Physics.h"
+
+/*
+
+
+ */
+
+class BorisPusher {
+    
+public:
+    BorisPusher(const PartData &ref);
+    BorisPusher();
+    void initialise(const PartData *ref);
+    
+    void kick(const Vector_t &R, Vector_t &P,
+              const Vector_t &Ef, const Vector_t &Bf,
+              const double &dt) const;
+    
+    void kick(const Vector_t &R, Vector_t &P,
+              const Vector_t &Ef, const Vector_t &Bf,
+              const double &dt, const double &mass,
+              const double &charge) const;
+    
+    void push(Vector_t &R, const Vector_t &P, const double &dt) const;
+
+private:
+    const PartData *itsReference;
+};
+
+inline BorisPusher::BorisPusher(const PartData &ref):
+    itsReference(&ref)
+{ }
+
+inline BorisPusher::BorisPusher():
+    itsReference(NULL)
+{ }
+
+inline void BorisPusher::initialise(const PartData *ref)
+{
+    itsReference = ref;
+}
+
+inline void BorisPusher::kick(const Vector_t &R, Vector_t &P,
+                              const Vector_t &Ef, const Vector_t &Bf,
+                              const double &dt) const
+{
+    kick(R, P, Ef, Bf, dt, itsReference->getM(), itsReference->getQ());
+}
+
+
+inline void BorisPusher::kick(const Vector_t &R, Vector_t &P,
+                              const Vector_t &Ef, const Vector_t &Bf,
+                              const double &dt, const double &mass,
+                              const double &charge) const
+{
+    // 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
+
+    using Physics::c;
+
+    // Half step E
+    P += 0.5 * dt * charge * c / mass * Ef;
+
+    // Full step B
+
+    /*
+        LF
+    double const gamma = sqrt(1.0 + dot(P, P));
+    Vector_t const t = dt * charge * c * c / (gamma * mass) * Bf;
+    P +=  cross(P, t);
+    */
+
+    double const gamma = sqrt(1.0 + dot(P, P));
+    Vector_t const t = 0.5 * dt * charge * c * c / (gamma * mass) * Bf;
+    Vector_t const w = P + cross(P, t);
+    Vector_t const s = 2.0 / (1.0 + dot(t, t)) * t;
+    P += cross(w, s);
+
+    /* a poor Leap-Frog
+        P += 1.0 * dt * charge * c * c / (gamma * mass) * cross(P,Bf);
+    */
+
+    // Half step E
+    P += 0.5 * dt * charge * c / mass * Ef;
+}
+
+
+inline void BorisPusher::push(Vector_t &R, const Vector_t &P, const double &/* dt */) const {
+    /** \f[ \vec{x}_{n+1/2} = \vec{x}_{n} + \frac{1}{2}\vec{v}_{n-1/2}\quad (= \vec{x}_{n} + \frac{\Delta t}{2} \frac{\vec{\beta}_{n-1/2}\gamma_{n-1/2}}{\gamma_{n-1/2}}) \f]
+     *
+     * \code
+     * R[i] += 0.5 * P[i] * recpgamma;
+     * \endcode
+     */
+    R += 0.5 * P / sqrt(1.0 + dot(P, P));
+}
+
+#endif
\ No newline at end of file

src/Steppers/CMakeLists.txt

+set (_SRCS
+)
+
+include_directories (
+  ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+add_opal_sources(${_SRCS})
+
+set (HDRS
+    BorisPusher.h
+    LF2.h
+    LF2.hpp
+    RK4.h
+    RK4.hpp
+    Stepper.h
+    Steppers.h
+)
+
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Steppers")
\ No newline at end of file

src/Steppers/LF2.h

+#ifndef LF2_H
+#define LF2_H
+
+#include "Stepper.h"
+#include "Physics/Physics.h"
+
+/// Leap-Frog 2nd order
+template <typename FieldFunction, typename ... Arguments>
+class LF2 : public Stepper<FieldFunction, Arguments...> {
+    
+public:
+    
+    LF2(const FieldFunction& fieldfunc) : Stepper<FieldFunction, Arguments ...>(fieldfunc) { }
+    
+    bool advance(PartBunch* bunch,
+                 const size_t& i,
+                 const double& t,
+                 const double dt,
+                 Arguments& ... args) const;
+    
+private:
+    
+    void push_m(Vector_t& R, const Vector_t& P, const double& h) const;
+    
+    bool kick_m(PartBunch* bunch, const size_t& i,
+                const double& t, const double& h,
+                Arguments& ... args) const;
+};
+
+#include "LF2.hpp"
+
+#endif

src/Steppers/LF2.hpp

+template <typename FieldFunction, typename ... Arguments>
+bool LF2<FieldFunction, Arguments ...>::advance(PartBunch* bunch,
+                                                const size_t& i,
+                                                const double& t,
+                                                const double dt,
+                                                Arguments& ... args) const
+{
+    bool flagNoDeletion = true;
+    
+    // push for first LF2 half step
+    push_m(bunch->R[i], bunch->P[i], 0.5 * dt * 1.0e-9);  // ns --> s
+    
+    flagNoDeletion = kick_m(bunch, i, t, dt * 1.0e-9, args ...);
+    
+    // push for second LF2 half step
+    push_m(bunch->R[i], bunch->P[i], 0.5 * dt * 1.0e-9);  // ns --> s
+    
+    return flagNoDeletion;
+}
+
+
+template <typename FieldFunction, typename ... Arguments>
+void LF2<FieldFunction, Arguments ...>::push_m(Vector_t& R, const Vector_t& P,
+                                               const double& h) const
+{
+    double const gamma = sqrt(1.0 + dot(P, P));
+    double const c_gamma = Physics::c / gamma;
+    Vector_t const v = P * c_gamma;
+    R += h * v;
+    
+    
+//     bunch->setT(bunch->getT() + h);
+// 
+//     // Path length update
+//     double const dotP = dot(bunch->P[0], bunch->P[0]);
+//     double const gamma = sqrt(1.0 + dotP);
+//     PathLength_m += h * sqrt(dotP) * Physics::c / gamma;
+//     bunch->setLPath(PathLength_m);
+}
+
+
+template <typename FieldFunction, typename ... Arguments>
+bool LF2<FieldFunction, Arguments ...>::kick_m(PartBunch* bunch, const size_t& i,
+                                               const double& t, const double& h,
+                                               Arguments& ... args) const
+{
+    Vector_t externalE = Vector_t(0.0, 0.0, 0.0);
+    Vector_t externalB = Vector_t(0.0, 0.0, 0.0);
+    
+    bool outOfBound = this->fieldfunc_m(t, i, externalE, externalB, args ...);
+    
+    if ( outOfBound )
+        return false;
+    
+    
+    double const q = bunch->Q[0] / Physics::q_e; // For now all particles have the same charge
+    double const M = bunch->M[0] * 1.0e9; // For now all particles have the same rest energy
+    double const h12Halfqc_M = 0.5 * h * q * Physics::c / M;
+    double const h12Halfqcc_M = h12Halfqc_M * Physics::c;
+    
+    // Half step E
+    bunch->P[i] += h12Halfqc_M * bunch->Ef[i];
+
+    // Full step B
+    double const gamma = sqrt(1.0 + dot(bunch->P[i], bunch->P[i]));
+    Vector_t const r = h12Halfqcc_M * bunch->Bf[i] / gamma;
+    Vector_t const w = bunch->P[i] + cross(bunch->P[i], r);
+    Vector_t const s = 2.0 / (1.0 + dot(r, r)) * r;
+    bunch->P[i] += cross(w, s);
+
+    // Half step E
+    bunch->P[i] += h12Halfqc_M * bunch->Ef[i];
+    
+    return true;
+}

src/Steppers/RK4.h

+#ifndef RK4_H
+#define RK4_H
+
+#include "Stepper.h"
+#include "Physics/Physics.h"
+
+/// 4-th order Runnge-Kutta stepper
+template <typename FieldFunction, typename ... Arguments>
+class RK4 : public Stepper<FieldFunction, Arguments...> {
+    
+public:
+    
+    RK4(const FieldFunction& fieldfunc) : Stepper<FieldFunction, Arguments ...>(fieldfunc) { }
+
+    bool advance(PartBunch* bunch,
+                 const size_t& i,
+                 const double& t,
+                 const double dt,
+                 Arguments& ... args) const;
+    
+private:
+    /**
+     * 
+     *
+     * @param y
+     * @param t
+     * @param yp
+     * @param Pindex
+     *
+     * @return
+     */
+    bool derivate_m(PartBunch* bunch,
+                    double *y,
+                    const double& t,
+                    double* yp,
+                    const size_t& i,
+                    Arguments& ... args) const;
+    
+    
+    void copyTo(const Vector_t& R, const Vector_t& P, double* x) const;
+    
+    void copyFrom(Vector_t& R, Vector_t& P, double* x) const;
+    
+    const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
+    const double c_mtns = Physics::c * 1.0e-9;  // m/s --> m/ns
+    const double mass_coeff = 1.0e18 * Physics::q_e / Physics::c / Physics::c; // from GeV/c^2 to basic unit: GV*C*s^2/m^2
+};
+
+#include "RK4.hpp"
+
+#endif