src/AbsBeamline/PluginElement.cpp

@@ -19,7 +19,7 @@
 #include "AbsBeamline/PluginElement.h"
 
 #include "AbsBeamline/BeamlineVisitor.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Physics/Physics.h"
 #include "Physics/Units.h"
 #include "Structure/LossDataSink.h"
@@ -77,6 +77,14 @@ bool PluginElement::apply(
     return false;
 }
 
+bool PluginElement::apply(
+        const Vector_t<double, 3>& R, const Vector_t<double, 3>& P, const double& t,
+        Vector_t<double, 3>& E, Vector_t<double, 3>& B) {
+    return false;
+}
+
+
+
 bool PluginElement::applyToReferenceParticle(
     const Vector_t<double, 3>&, const Vector_t<double, 3>&, const double&, Vector_t<double, 3>&,
     Vector_t<double, 3>&) {

src/AbsBeamline/PluginElement.h

@@ -47,6 +47,11 @@ public:
     virtual bool apply(
         const size_t& i, const double& t, Vector_t<double, 3>& E, Vector_t<double, 3>& B) override;
 
+    virtual bool apply(
+        const Vector_t<double, 3>& R, const Vector_t<double, 3>& P, const double& t,
+        Vector_t<double, 3>& E, Vector_t<double, 3>& B);
+
+    
     virtual bool applyToReferenceParticle(
         const Vector_t<double, 3>& R, const Vector_t<double, 3>& P, const double& t,
         Vector_t<double, 3>& E, Vector_t<double, 3>& B) override;

src/AbsBeamline/Probe.cpp

@@ -17,7 +17,7 @@
 //
 #include "AbsBeamline/Probe.h"
 #include "AbsBeamline/BeamlineVisitor.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Physics/Physics.h"
 #include "Physics/Units.h"
 #include "Structure/LossDataSink.h"

src/AbsBeamline/RFCavity.cpp

@@ -21,7 +21,7 @@
 #include <boost/filesystem.hpp>
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "Fields/Fieldmap.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Physics/Units.h"
 #include "Steppers/BorisPusher.h"
 #include "Utilities/GeneralClassicException.h"

src/AbsBeamline/Ring.cpp

@@ -35,7 +35,7 @@
 
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "Fields/EMField.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Physics/Physics.h"
 #include "Structure/LossDataSink.h"
 

src/AbsBeamline/ScalingFFAMagnet.cpp

@@ -26,10 +26,9 @@
  */
 
 #include <cmath>
-
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "AbsBeamline/ScalingFFAMagnet.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Physics/Units.h"
 ScalingFFAMagnet::ScalingFFAMagnet(const std::string& name)
     : Component(name), planarArcGeometry_m(1., 1.), dummy(), endField_m(nullptr) {

src/AbsBeamline/ScalingFFAMagnet.h

@@ -30,6 +30,29 @@
 #include "AbsBeamline/EndFieldModel/EndFieldModel.h"
 #include "AbsBeamline/Component.h"
 
+#ifdef __CUDACC__
+#pragma push_macro("__cpp_consteval")
+#pragma push_macro("_NODISCARD")
+#pragma push_macro("__builtin_LINE")
+
+#define __cpp_consteval 201811L
+
+#ifdef _NODISCARD
+    #undef _NODISCARD
+    #define _NODISCARD
+#endif
+
+#define consteval constexpr
+
+#include <source_location>
+
+#undef consteval
+#pragma pop_macro("__cpp_consteval")
+#pragma pop_macro("_NODISCARD")
+#else
+#include <source_location>
+#endif
+
 #ifndef ABSBEAMLINE_ScalingFFAMagnet_H
 #define ABSBEAMLINE_ScalingFFAMagnet_H
 

src/AbsBeamline/Solenoid.cpp

@@ -21,7 +21,7 @@
 #include "AbsBeamline/Solenoid.h"
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "Fields/Fieldmap.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Physics/Physics.h"
 
 #include <fstream>

src/AbsBeamline/Solenoid.h

@@ -21,6 +21,29 @@
 //
 // ------------------------------------------------------------------------
 
+#ifdef __CUDACC__
+#pragma push_macro("__cpp_consteval")
+#pragma push_macro("_NODISCARD")
+#pragma push_macro("__builtin_LINE")
+
+#define __cpp_consteval 201811L
+
+#ifdef _NODISCARD
+    #undef _NODISCARD
+    #define _NODISCARD
+#endif
+
+#define consteval constexpr
+
+#include <source_location>
+
+#undef consteval
+#pragma pop_macro("__cpp_consteval")
+#pragma pop_macro("_NODISCARD")
+#else
+#include <source_location>
+#endif
+
 #include "AbsBeamline/Component.h"
 #include "Algorithms/CoordinateSystemTrafo.h"
 

src/AbsBeamline/TravelingWave.cpp

@@ -18,7 +18,7 @@
 #include "AbsBeamline/TravelingWave.h"
 
 #include "AbsBeamline/BeamlineVisitor.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Fields/Fieldmap.h"
 #include "Physics/Units.h"
 
@@ -235,6 +235,13 @@ void TravelingWave::initialise(PartBunch_t* bunch, double& startField, double& e
         - Physics::two_pi * ((numCells_m - 1) * mode_m - std::floor((numCells_m - 1) * mode_m));
 }
 
+void TravelingWave::initialise(PartBunch_t* bunch, std::shared_ptr<AbstractTimeDependence> freq_atd,
+                               std::shared_ptr<AbstractTimeDependence> ampl_atd,
+                               std::shared_ptr<AbstractTimeDependence> phase_atd) {
+    
+}
+
+
 void TravelingWave::finalise() {
 }
 

src/AbsBeamline/TravelingWave.h

@@ -68,6 +68,10 @@ public:
 
     virtual void initialise(PartBunch_t* bunch, double& startField, double& endField) override;
 
+    virtual void initialise(PartBunch_t* bunch, std::shared_ptr<AbstractTimeDependence> freq_atd,
+        std::shared_ptr<AbstractTimeDependence> ampl_atd,
+        std::shared_ptr<AbstractTimeDependence> phase_atd) override;
+    
     virtual void finalise() override;
 
     virtual bool bends() const override;

src/AbsBeamline/VerticalFFAMagnet.h

@@ -10,7 +10,7 @@
 #include "AbsBeamline/Component.h"
 #include "BeamlineGeometry/StraightGeometry.h"
 #include "Fields/BMultipoleField.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 
 #ifndef ABSBEAMLINE_VerticalFFAMagnet_H
 #define ABSBEAMLINE_VerticalFFAMagnet_H

src/AbstractObjects/Object.h

@@ -25,6 +25,29 @@
 #include <string>
 #include <vector>
 
+#ifdef __CUDACC__
+#pragma push_macro("__cpp_consteval")
+#pragma push_macro("_NODISCARD")
+#pragma push_macro("__builtin_LINE")
+
+#define __cpp_consteval 201811L
+
+#ifdef _NODISCARD
+    #undef _NODISCARD
+    #define _NODISCARD
+#endif
+
+#define consteval constexpr
+
+#include <source_location>
+
+#undef consteval
+#pragma pop_macro("__cpp_consteval")
+#pragma pop_macro("_NODISCARD")
+#else
+#include <source_location>
+#endif
+
 class Invalidator;
 class Parser;
 class Statement;
@@ -117,7 +140,7 @@ public:
 
     /// Print the object.
     //  Print a OPAL-readable image of [b]this[/b] on the given output stream.
-    virtual void print(std::ostream &) const;
+    virtual void print(std::ostream  &) const;
 
     virtual void printValue(std::ostream &) const;
 

src/AbstractObjects/OpalData.cpp

@@ -28,7 +28,7 @@
 #include "AbstractObjects/ObjectFunction.h"
 #include "AbstractObjects/Table.h"
 #include "AbstractObjects/ValueDefinition.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Attributes/Attributes.h"
 #include "OpalParser/FileStream.h"
 #include "OpalParser/OpalParser.h"
@@ -786,4 +786,4 @@ void OpalData::storeArguments(int argc, char* argv[]) {
 
 std::vector<std::string> OpalData::getArguments() {
     return p->arguments_m;
-}
\ No newline at end of file
+}

src/AbstractObjects/OpalData.h

@@ -22,6 +22,29 @@
 #ifndef OPAL_OpalData_HH
 #define OPAL_OpalData_HH
 
+#ifdef __CUDACC__
+#pragma push_macro("__cpp_consteval")
+#pragma push_macro("_NODISCARD")
+#pragma push_macro("__builtin_LINE")
+
+#define __cpp_consteval 201811L
+
+#ifdef _NODISCARD
+    #undef _NODISCARD
+    #define _NODISCARD
+#endif
+
+#define consteval constexpr
+
+#include <source_location>
+
+#undef consteval
+#pragma pop_macro("__cpp_consteval")
+#pragma pop_macro("_NODISCARD")
+#else
+#include <source_location>
+#endif
+
 #include <iosfwd>
 #include <map>
 #include <stack>

src/Algorithms/DistributionMoments.cpp

@@ -44,20 +44,30 @@ DistributionMoments::DistributionMoments() {
     notCentMoments_m.resize(6, 6, false);
 }
 
-
-
 void DistributionMoments::computeMeans(ippl::ParticleAttrib<Vector_t<double,3>>::view_type&  Rview,
                                          ippl::ParticleAttrib<Vector_t<double,3>>::view_type&  Pview,
                                          ippl::ParticleAttrib<double>::view_type&  Mview,
-                                         size_t Np) {
-    const int Dim = 3;
+                                         size_t Np,
+                                         size_t Nlocal) {
+    /*
+     Np is the total number of particles (reduced over ranks). In this function, it is only used to 
+     average over the number of total particles. For an empty simulation, this leads to divisions by
+     0. Since, however, some of the computed moments might be needed regardless, we compute them but
+     need to make sure that we do not divide by 0. 
+     Solution: Set Np to 1 if it is 0.
+     */
+    Np = (Np == 0) ? 1 : Np;
 
+    const int Dim = 3;
     double loc_centroid[2 * Dim]        = {};
     double centroid[2 * Dim]        = {};
     double loc_Ekin, loc_gamma, loc_gammaz, gammaz;
 
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
     Kokkos::parallel_reduce(
-                                    "calc moments of particle distr.", ippl::getRangePolicy(Rview),
+                                    "calc moments of particle distr.", Nlocal,
                 KOKKOS_LAMBDA(
                     const int k, double& ekin, double& gamma, double& gammaz) {
                     double gamma0 = 0.0;
@@ -78,7 +88,7 @@ void DistributionMoments::computeMeans(ippl::ParticleAttrib<Vector_t<double,3>>:
 
     for (unsigned i = 0; i < 2 * Dim; ++i) {
             Kokkos::parallel_reduce(
-                                    "calc moments of particle distr.", ippl::getRangePolicy(Rview),
+                                    "calc moments of particle distr.", Nlocal,
                 KOKKOS_LAMBDA(
                     const int k, double& cent) {
                     double part[2 * Dim];
@@ -126,9 +136,12 @@ void DistributionMoments::computeMeans(ippl::ParticleAttrib<Vector_t<double,3>>:
 void DistributionMoments::computeMoments(ippl::ParticleAttrib<Vector_t<double,3>>::view_type&  Rview,
                                          ippl::ParticleAttrib<Vector_t<double,3>>::view_type&  Pview,
                                          ippl::ParticleAttrib<double>::view_type&  Mview,
-                                         size_t Np) {
+                                         size_t Np,
+                                         size_t Nlocal) {
+    Np = (Np == 0) ? 1 : Np; // Explanation: see DistributionMoments::computeMeans implementation
+
     reset();
-    computeMeans(Rview, Pview, Mview, Np);
+    computeMeans(Rview, Pview, Mview, Np, Nlocal);
 
     double meanR_loc[Dim]       = {};
     double meanP_loc[Dim]       = {};
@@ -146,7 +159,7 @@ void DistributionMoments::computeMoments(ippl::ParticleAttrib<Vector_t<double,3>
 
     for (unsigned i = 0; i < 2 * Dim; ++i) {
             Kokkos::parallel_reduce(
-                                    "calc moments of particle distr.", ippl::getRangePolicy(Rview),
+                                    "calc moments of particle distr.", Nlocal,
                 KOKKOS_LAMBDA(
                     const int k, double& mom0, double& mom1, double& mom2,
                     double& mom3, double& mom4, double& mom5) {
@@ -171,14 +184,13 @@ void DistributionMoments::computeMoments(ippl::ParticleAttrib<Vector_t<double,3>
                 Kokkos::Sum<double>(loc_moment[i][5]));
             Kokkos::fence();
      }
-    ippl::Comm->barrier();
 
     MPI_Allreduce(
             loc_moment, moment, 2 * Dim * 2 * Dim, MPI_DOUBLE, MPI_SUM, ippl::Comm->getCommunicator());
 
     for (unsigned i = 0; i < 2 * Dim; i++) {
             for (unsigned j = 0; j < 2 * Dim; j++) {
-                moments_m(i,j) = moment[i][j] / (Np-1);
+                moments_m(i,j) = moment[i][j] / Np;
             }
      }
 
@@ -193,7 +205,7 @@ void DistributionMoments::computeMoments(ippl::ParticleAttrib<Vector_t<double,3>
    // compute non-central moments
    for (unsigned i = 0; i < 2 * Dim; ++i) {
             Kokkos::parallel_reduce(
-                                    "calc moments of particle distr.", ippl::getRangePolicy(Rview),
+                                    "calc moments of particle distr.", Nlocal,
                 KOKKOS_LAMBDA(
                     const int k, double& mom0, double& mom1, double& mom2,
                     double& mom3, double& mom4, double& mom5) {
@@ -221,7 +233,7 @@ void DistributionMoments::computeMoments(ippl::ParticleAttrib<Vector_t<double,3>
     ippl::Comm->barrier();
 
     Kokkos::parallel_reduce(
-                "calc moments of particle distr.", ippl::getRangePolicy(Pview),
+                "calc moments of particle distr.", Nlocal,
                 KOKKOS_LAMBDA(
                     const int k, double& ekin) {
                     double gamma0 = 0;
@@ -279,11 +291,12 @@ void DistributionMoments::computeMoments(ippl::ParticleAttrib<Vector_t<double,3>
     }
 
     double betaGamma = std::sqrt(std::pow(meanGamma_m, 2) - 1.0);
-    geometricEps_m = normalizedEps_m / Vector_t(betaGamma);
+    geometricEps_m = normalizedEps_m / Vector_t<double,3>(betaGamma);
 
 }
 
-void DistributionMoments::computeMinMaxPosition(ippl::ParticleAttrib<Vector_t<double,3>>::view_type& Rview){
+void DistributionMoments::computeMinMaxPosition(ippl::ParticleAttrib<Vector_t<double,3>>::view_type& Rview, size_t Nlocal)
+{
     const int Dim = 3;
 
     double rmax_loc[Dim];
@@ -291,9 +304,15 @@ void DistributionMoments::computeMinMaxPosition(ippl::ParticleAttrib<Vector_t<do
     double rmax[Dim];
     double rmin[Dim];
 
+    for(int i=0; i<Dim; i++){
+        rmin_loc[i] = 0.;
+        rmax_loc[i] = 0.;
+    }
+
     for (unsigned d = 0; d < Dim; ++d) {
+        if (Nlocal > 0) {
             Kokkos::parallel_reduce(
-                "rel max", ippl::getRangePolicy(Rview),
+                "rel max", Nlocal,
                 KOKKOS_LAMBDA(const int i, double& mm) {
                     double tmp_vel = Rview(i)[d];
                     mm             = tmp_vel > mm ? tmp_vel : mm;
@@ -307,6 +326,7 @@ void DistributionMoments::computeMinMaxPosition(ippl::ParticleAttrib<Vector_t<do
                     mm             = tmp_vel < mm ? tmp_vel : mm;
                 },
                 Kokkos::Min<double>(rmin_loc[d]));
+         }
      }
      Kokkos::fence();
      MPI_Allreduce(rmax_loc, rmax, Dim, MPI_DOUBLE, MPI_MAX, ippl::Comm->getCommunicator());
@@ -724,7 +744,10 @@ void DistributionMoments::computeMeanKineticEnergy() {
 void DistributionMoments::computeDebyeLength(ippl::ParticleAttrib<Vector_t<double,3>>::view_type&  Rview,
                                          ippl::ParticleAttrib<Vector_t<double,3>>::view_type&  Pview,
                                          size_t Np,
+                                         size_t Nlocal,
                                          double density){
+    Np = (Np == 0) ? 1 : Np; // Explanation: see DistributionMoments::computeMeans implementation
+
     resetPlasmaParameters();
     double loc_avgVel[3] = {};
     double avgVel[3] = {};
@@ -732,7 +755,7 @@ void DistributionMoments::computeDebyeLength(ippl::ParticleAttrib<Vector_t<doubl
 
     // From P in \beta\gamma to get v in m/s: v = (P*c)/\gamma
     Kokkos::parallel_reduce(
-                "calc moments of particle distr.", ippl::getRangePolicy(Rview),
+                "calc moments of particle distr.", Nlocal,
                 KOKKOS_LAMBDA(
                     const int k, double& mom0, double& mom1, double& mom2){
 
@@ -773,7 +796,7 @@ void DistributionMoments::computeDebyeLength(ippl::ParticleAttrib<Vector_t<doubl
     }*/
 
     Kokkos::parallel_reduce(
-                "calc moments of particle distr.", ippl::getRangePolicy(Rview),
+                "calc moments of particle distr.", Nlocal,
                 KOKKOS_LAMBDA(
                     const int k, double& mom0){
 

src/Algorithms/DistributionMoments.h

@@ -18,6 +18,29 @@
 #ifndef DISTRIBUTIONMOMENTS_H
 #define DISTRIBUTIONMOMENTS_H
 
+#ifdef __CUDACC__
+#pragma push_macro("__cpp_consteval")
+#pragma push_macro("_NODISCARD")
+#pragma push_macro("__builtin_LINE")
+
+#define __cpp_consteval 201811L
+
+#ifdef _NODISCARD
+    #undef _NODISCARD
+    #define _NODISCARD
+#endif
+
+#define consteval constexpr
+
+#include <source_location>
+
+#undef consteval
+#pragma pop_macro("__cpp_consteval")
+#pragma pop_macro("_NODISCARD")
+#else
+#include <source_location>
+#endif
+
 #include "Ippl.h"
 #include <Kokkos_Core.hpp>
 #include "Algorithms/BoostMatrix.h"
@@ -45,12 +68,14 @@ public:
     void computeMoments(ippl::ParticleAttrib<Vector_t<double,3>>::view_type& Rview,
                         ippl::ParticleAttrib<Vector_t<double,3>>::view_type& Pview,
                         ippl::ParticleAttrib<double>::view_type& Mview,
-                        size_t Np);
-    void computeMinMaxPosition(ippl::ParticleAttrib<Vector_t<double,3>>::view_type& Rview);
+                        size_t Np,
+                        size_t Nlocal);
+    void computeMinMaxPosition(ippl::ParticleAttrib<Vector_t<double,3>>::view_type& Rview, size_t Nlcoal);
     void computeMeanKineticEnergy();
     void computeDebyeLength(ippl::ParticleAttrib<Vector_t<double,3>>::view_type& Rview,
                         ippl::ParticleAttrib<Vector_t<double,3>>::view_type& Pview,
                         size_t Np,
+                        size_t Nlocal,
                         double density);
     void computePlasmaParameter(double);
 
@@ -95,12 +120,13 @@ public:
     double getTotalMass() const;
     double getTotalNumParticles() const;
 
-private:
-    bool isParticleExcluded(const OpalParticle&) const;
     void computeMeans(ippl::ParticleAttrib<Vector_t<double,3>>::view_type&  Rview,
                                          ippl::ParticleAttrib<Vector_t<double,3>>::view_type&  Pview,
                                          ippl::ParticleAttrib<double>::view_type&  Mview,
-                                         size_t Np);
+                                         size_t Np,
+                                         size_t Nlocal);
+private:
+    bool isParticleExcluded(const OpalParticle&) const;
 
     //template <class InputIt>
     //void computeMeans(const InputIt&, const InputIt&);

src/Algorithms/OrbitThreader.cpp

@@ -134,14 +134,13 @@ void OrbitThreader::execute() {
     errorFlag_m = EVERYTHINGFINE;
 
     *gmsg << "OrbitThreader dt_m= " << dt_m << endl;
-
-
+    
     do {
         checkElementLengths(elementSet);
         if (containsCavity(elementSet)) {
             autophaseCavities(elementSet, visitedElements);
         }
-
+        
         double initialS              = pathLength_m;
         Vector_t<double, 3> initialR = r_m;
         Vector_t<double, 3> initialP = p_m;
@@ -149,6 +148,9 @@ void OrbitThreader::execute() {
 
         integrate(elementSet, maxDistance);
 
+        *gmsg << "OrbitThreader maxDistance= " << maxDistance << endl;
+        *gmsg << "OrbitThreader #elements  = " << elementSet.size() << endl;
+        
         registerElement(elementSet, initialS, initialR, initialP);
 
         if (errorFlag_m == HITMATERIAL) {
@@ -251,6 +253,7 @@ void OrbitThreader::integrate(const IndexMap::value_t& activeSet, double /*maxDr
         const Vector<double, 3> d = r_m - oldR;
 
         pathLength_m += std::copysign(euclidean_norm(d), dt_m);
+
         ++currentStep_m;
         time_m += dt_m;
 
@@ -439,7 +442,6 @@ void OrbitThreader::computeBoundingBox() {
     FieldList allElements         = itsOpalBeamline_m.getElementByType(ElementType::ANY);
     FieldList::iterator it        = allElements.begin();
     const FieldList::iterator end = allElements.end();
-
     for (; it != end; ++it) {
         if (it->getElement()->getType() == ElementType::MARKER) {
             continue;
@@ -452,26 +454,25 @@ void OrbitThreader::computeBoundingBox() {
 
 void OrbitThreader::updateBoundingBoxWithCurrentPosition() {
     Vector_t<double, 3> dR = Physics::c * dt_m * p_m / Util::getGamma(p_m);
-
-    /// \todo needs to be fixed
-    /*
-      for (const Vector_t<double, 3>& pos : {r_m - 10 * dR, r_m + 10 * dR}) {
+    std::array<Vector_t<double, 3>, 2> positions = {r_m - 10 * dR, r_m + 10 * dR};
+   
+    for (const Vector_t<double, 3>& pos : positions) {
         globalBoundingBox_m.enlargeToContainPosition(pos);
     }
-    */
+    
 }
 
 double OrbitThreader::computeDriftLengthToBoundingBox(
     const std::set<std::shared_ptr<Component>>& elements, const Vector_t<double, 3>& position,
     const Vector_t<double, 3>& direction) const {
-    /*
+
     if (elements.empty()
         || (elements.size() == 1 && (*elements.begin())->getType() == ElementType::DRIFT)) {
         boost::optional<Vector_t<double, 3>> intersectionPoint =
             globalBoundingBox_m.getIntersectionPoint(position, direction);
-
-        return intersectionPoint ? euclidean_norm(intersectionPoint.get() - position) : 10.0;
+        const Vector_t<double, 3> r = intersectionPoint.get() - position;
+        return intersectionPoint ? euclidean_norm(r) : 10.0;
     }
-    */
+
     return std::numeric_limits<double>::max();
 }
\ No newline at end of file

src/Algorithms/ParallelTracker.cpp

@@ -71,14 +71,15 @@ ParallelTracker::ParallelTracker(
       zstart_m(0.0),
       dtCurrentTrack_m(0.0),
       minStepforReBin_m(-1),
-      repartFreq_m(-1),
+      repartFreq_m(0),
       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")),
       PluginElemTimer_m(IpplTimings::getTimer("PluginElements")),
-      BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")) {
+      BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
+      OrbThreader_m(IpplTimings::getTimer("OrbThreader")) {
 }
 
 ParallelTracker::ParallelTracker(
@@ -96,22 +97,21 @@ ParallelTracker::ParallelTracker(
       zstart_m(zstart),
       dtCurrentTrack_m(0.0),
       minStepforReBin_m(-1),
-      repartFreq_m(-1),
+      repartFreq_m(0),
       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")) {
-
-    *gmsg << "* ParallelTracker zstop.size()= " << zstop.size() << endl;
-
-    for (unsigned int i = 0; i < zstop.size(); ++i) {
-        stepSizes_m.push_back(dt[i], zstop[i], maxSteps[i]);
-    }
+      BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
+      OrbThreader_m(IpplTimings::getTimer("OrbThreader")) {
+    
+      for (unsigned int i = 0; i < zstop.size(); ++i) {
+          stepSizes_m.push_back(dt[i], zstop[i], maxSteps[i]);
+      }
 
-    stepSizes_m.sortAscendingZStop();
-    stepSizes_m.resetIterator();
+      stepSizes_m.sortAscendingZStop();
+      stepSizes_m.resetIterator();
 }
 
 ParallelTracker::~ParallelTracker() {
@@ -348,13 +348,23 @@ void ParallelTracker::execute() {
     *gmsg << "itsBunch_m->RefPartR_m= " << itsBunch_m->RefPartR_m << endl;
     *gmsg << "itsBunch_m->RefPartP_m= " << itsBunch_m->RefPartP_m << endl;
 
+    IpplTimings::startTimer(OrbThreader_m);
+
+    bool const psDump0 = 0;
+    bool const statDump0 = 0;
+
+    writePhaseSpace(0, psDump0, statDump0);
+    msg << level2 << "Dump initial phase space" << endl;
+
+
     OrbitThreader oth(
         itsReference, itsBunch_m->RefPartR_m, itsBunch_m->RefPartP_m, pathLength_m, -rmin(2),
         itsBunch_m->getT(), (back_track ? -minTimeStep : minTimeStep), stepSizes_m,
         itsOpalBeamline_m);
 
     oth.execute();
-
+    IpplTimings::stopTimer(OrbThreader_m);
+    
     BoundingBox globalBoundingBox = oth.getBoundingBox();
 
     numParticlesInSimulation_m = itsBunch_m->getTotalNum();
@@ -368,7 +378,6 @@ void ParallelTracker::execute() {
     OPALTimer::Timer myt1;
     *gmsg << "* Track start at: " << myt1.time() << ", t= " << Util::getTimeString(time) << "; "
           << "zstart at: " << Util::getLengthString(pathLength_m) << endl;
-
     *gmsg << "* Executing ParallelTracker\n"
           << "* Initial dt = " << Util::getTimeString(itsBunch_m->getdT()) << "\n"
           << "* Max integration steps = " << stepSizes_m.getMaxSteps() << ", next step = " << step
@@ -383,8 +392,9 @@ void ParallelTracker::execute() {
     OpalData::getInstance()->setInPrepState(false);
 
     stepSizes_m.printDirect(*gmsg);
-
+    
     while (!stepSizes_m.reachedEnd()) {
+
         unsigned long long trackSteps = stepSizes_m.getNumSteps() + step;
         dtCurrentTrack_m              = stepSizes_m.getdT();
         changeDT(back_track);
@@ -394,23 +404,22 @@ void ParallelTracker::execute() {
             if (itsBunch_m->getTotalNum() > 0) {
                 itsBunch_m->get_bounds(rmin, rmax);
             }
-
+            // ADA
             timeIntegration1(pusher);
-
+            
             computeSpaceChargeFields(step);
             
             // \todo for a drift we can neglect that 
             // computeExternalFields(oth);
 
             timeIntegration2(pusher);
-
+            
             selectDT(back_track);
             // \todo emitParticles(step);
             //selectDT(back_track);
 
             itsBunch_m->incrementT();
 
-            
             if (itsBunch_m->getT() > 0.0 || itsBunch_m->getdT() < 0.0) {
                 updateReference(pusher);
             }
@@ -421,7 +430,7 @@ void ParallelTracker::execute() {
             }
 
             itsBunch_m->set_sPos(pathLength_m);
-            
+
             // if (hasEndOfLineReached(globalBoundingBox)) break;
   
             bool const psDump =
@@ -459,7 +468,7 @@ void ParallelTracker::execute() {
 
     writePhaseSpace((step + 1), psDump, statDump);
 
-    msg << level2 << "Dump phase space of last step" << endl;
+    *gmsg << "* Dump phase space of last step" << endl;
 
     itsOpalBeamline_m.switchElementsOff();
 
@@ -508,6 +517,7 @@ void ParallelTracker::timeIntegration2(BorisPusher& pusher) {
     // switchElements();
     pushParticles(pusher);
 
+    
     auto dtview  = itsBunch_m->getParticleContainer()->dt.getView();
     double newdT = itsBunch_m->getdT();
 
@@ -516,7 +526,7 @@ void ParallelTracker::timeIntegration2(BorisPusher& pusher) {
                          KOKKOS_LAMBDA(const int i) {
                              dtview(i) = newdT;
                          });                     
-
+    
     IpplTimings::stopTimer(timeIntegrationTimer2_m);
 }
 
@@ -572,6 +582,7 @@ void ParallelTracker::computeSpaceChargeFields(unsigned long long step) {
 
      */
 
+    
     const matrix_t                rot = referenceToBeamCSTrafo.getRotationMatrix();
     const ippl::Vector<double, 3> org = referenceToBeamCSTrafo.getOrigin();
 
@@ -627,7 +638,7 @@ void ParallelTracker::computeSpaceChargeFields(unsigned long long step) {
             prod_boost_vector(boost::numeric::ublas::trans(rotationMatrix_m)
 
      */
-
+    
     Kokkos::parallel_for(
                          "CSTrafo:transformTo", ippl::getRangePolicy(Rview),
                          KOKKOS_LAMBDA(const int i) {                           
@@ -753,18 +764,19 @@ void ParallelTracker::doBinaryRepartition() {
 
 void ParallelTracker::dumpStats(long long step, bool psDump, bool statDump) {
     OPALTimer::Timer myt2;
-    Inform msg("ParallelTracker ", *gmsg);
 
+    /*
     if (itsBunch_m->getGlobalTrackStep() % 1000 + 1 == 1000) {
-        msg << level1;
+        *gmsg << level1;
     } else if (itsBunch_m->getGlobalTrackStep() % 100 + 1 == 100) {
-        msg << level2;
+        *gmsg << level2;
     } else {
-        msg << level3;
+        *gmsg << level3;
     }
-
+    */
+    
     if (numParticlesInSimulation_m == 0) {
-        msg << myt2.time() << " "
+        *gmsg << "* " << myt2.time() << " "
             << "Step " << std::setw(6) << itsBunch_m->getGlobalTrackStep() << "; "
             << "   -- no emission yet --     "
             << "t= " << Util::getTimeString(itsBunch_m->getT()) << endl;
@@ -778,7 +790,7 @@ void ParallelTracker::dumpStats(long long step, bool psDump, bool statDump) {
             "ParallelTracker::dumpStats()",
             "there seems to be something wrong with the position of the bunch!");
     } else {
-        msg << myt2.time() << " "
+        *gmsg << "* " << myt2.time() << " "
             << "Step " << std::setw(6) << itsBunch_m->getGlobalTrackStep() << " "
             << "at " << Util::getLengthString(pathLength_m) << ", "
             << "t= " << Util::getTimeString(itsBunch_m->getT()) << ", "
@@ -790,7 +802,6 @@ void ParallelTracker::dumpStats(long long step, bool psDump, bool statDump) {
 
 void ParallelTracker::setOptionalVariables() {
 
-    Inform msg("ParallelTracker ", *gmsg);
     /*
     minStepforReBin_m = Options::minStepForRebin;
     RealVariable* br =
@@ -802,14 +813,14 @@ void ParallelTracker::setOptionalVariables() {
 
     // there is no point to do repartitioning with one node
     if (ippl::Comm->size() == 1) {
-        repartFreq_m = std::numeric_limits<unsigned int>::max();
+        repartFreq_m = std::numeric_limits<unsigned long long>::max();
     } else {
         repartFreq_m = Options::repartFreq * 100;
         RealVariable* rep =
             dynamic_cast<RealVariable*>(OpalData::getInstance()->find("REPARTFREQ"));
         if (rep)
             repartFreq_m = static_cast<int>(rep->getReal());
-        msg << level2 << "REPARTFREQ " << repartFreq_m << endl;
+        *gmsg  << "* REPARTFREQ " << repartFreq_m << endl;
     }
 }
 
@@ -833,8 +844,6 @@ void ParallelTracker::setTime() {
 }
 
 void ParallelTracker::writePhaseSpace(const long long /*step*/, bool psDump, bool statDump) {
-    extern Inform* gmsg;
-    Inform msg("OPAL ", *gmsg);
     Vector_t<double, 3> externalE, externalB;
     Vector_t<double, 3> FDext[2];  // FDext = {BHead, EHead, BRef, ERef, BTail, ETail}.
 
@@ -858,12 +867,14 @@ void ParallelTracker::writePhaseSpace(const long long /*step*/, bool psDump, boo
 
     if (statDump) {
         itsDataSink_m->dumpSDDS(itsBunch_m, FDext, -1.0);
-        msg << level3 << "* Wrote beam statistics." << endl;
+        *gmsg << "* Wrote beam statistics." << endl;
     }
 
-    /* \todo
-
     if (psDump && (itsBunch_m->getTotalNum() > 0)) {
+        // Write bunch to .h5 file.
+        itsDataSink_m->dumpH5(itsBunch_m, FDext);
+
+        /*
         // Write fields to .h5 file.
         const size_t localNum    = itsBunch_m->getLocalNum();
         double distToLastStop    = stepSizes_m.getFinalZStop() - pathLength_m;
@@ -910,9 +921,13 @@ void ParallelTracker::writePhaseSpace(const long long /*step*/, bool psDump, boo
         if (!statDump && !driftToCorrectPosition)
             itsBunch_m->calcBeamParameters();
 
-        msg << *itsBunch_m << endl;
+            msg << *itsBunch_m << endl;
+        
         itsDataSink_m->dumpH5(itsBunch_m, FDext);
-
+        */
+        
+        /*
+        
         if (driftToCorrectPosition) {
             if (localNum > 0) {
                 itsBunch_m->R = stashedR;
@@ -925,9 +940,10 @@ void ParallelTracker::writePhaseSpace(const long long /*step*/, bool psDump, boo
             itsBunch_m->calcBeamParameters();
         }
 
-        msg << level2 << "* Wrote beam phase space." << endl;
+        */
+        
+        *gmsg << level2 << "* Wrote beam phase space." << endl;
     }
-    */
 }
 
 void ParallelTracker::updateReference(const BorisPusher& pusher) {
@@ -977,7 +993,7 @@ void ParallelTracker::updateReferenceParticle(const BorisPusher& pusher) {
 void ParallelTracker::transformBunch(const CoordinateSystemTrafo& trafo) {
     const unsigned int localNum = itsBunch_m->getLocalNum();
     for (unsigned int i = 0; i < localNum; ++i) {
-        /*
+        /* \todo host device .... 
         itsBunch_m->R[i]  = trafo.transformTo(itsBunch_m->R[i]);
         itsBunch_m->P[i]  = trafo.rotateTo(itsBunch_m->P[i]);
         itsBunch_m->Ef[i] = trafo.rotateTo(itsBunch_m->Ef[i]);
@@ -987,11 +1003,21 @@ void ParallelTracker::transformBunch(const CoordinateSystemTrafo& trafo) {
 }
 
 void ParallelTracker::updateRefToLabCSTrafo() {
-    itsBunch_m->toLabTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
-    Vector_t<double, 3> R = itsBunch_m->RefPartR_m;
+    /*
+    Inform m("updateRefToLabCSTrafo ");
+    m << " initial RefPartR: " << itsBunch_m->RefPartR_m << endl;
+    m << " RefPartR after transformFrom( " << itsBunch_m->RefPartR_m << endl;
     itsBunch_m->toLabTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
     Vector_t<double, 3> P = itsBunch_m->RefPartP_m;
+    m << " CS " << itsBunch_m->toLabTrafo_m << endl;
+    m << " rotMat= " << itsBunch_m->toLabTrafo_m.getRotationMatrix() << endl;
+    m << " initial: path Lenght= " << pathLength_m << endl;
+    m << " dt= " << itsBunch_m->getdT() << " R=" << R << endl;
+    */
 
+    Vector_t<double, 3> R = itsBunch_m->toLabTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
+    Vector_t<double, 3> P = itsBunch_m->toLabTrafo_m.transformFrom(itsBunch_m->RefPartP_m);
+    
     pathLength_m += std::copysign(1, itsBunch_m->getdT()) * euclidean_norm(R);
 
     CoordinateSystemTrafo update(R, getQuaternion(P, Vector_t<double, 3>(0, 0, 1)));

src/Algorithms/ParallelTracker.h

@@ -95,7 +95,7 @@ class ParallelTracker : public Tracker {
     int minStepforReBin_m;
 
     // this variable controls the minimal number of steps until we repartition the particles
-    unsigned int repartFreq_m;
+    unsigned long long repartFreq_m;
 
     unsigned int emissionSteps_m;
 
@@ -107,7 +107,8 @@ class ParallelTracker : public Tracker {
     IpplTimings::TimerRef WakeFieldTimer_m;
     IpplTimings::TimerRef PluginElemTimer_m;
     IpplTimings::TimerRef BinRepartTimer_m;
-
+    IpplTimings::TimerRef OrbThreader_m;
+    
     std::set<ParticleMatterInteractionHandler*> activeParticleMatterInteractionHandlers_m;
     bool particleMatterStatus_m;
 
@@ -180,6 +181,18 @@ public:
     /// Apply the algorithm to a vertical FFA magnet.
     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();
@@ -221,8 +234,6 @@ private:
 
     /********************** END VARIABLES ***********************************/
 
-    void kickParticles(const BorisPusher& pusher);
-    void pushParticles(const BorisPusher& pusher);
     void updateReferenceParticle(const BorisPusher& pusher);
 
     void writePhaseSpace(const long long step, bool psDump, bool statDump);
@@ -237,22 +248,18 @@ private:
     void prepareSections();
 
     void timeIntegration1(BorisPusher& pusher);
-    void timeIntegration2(BorisPusher& pusher);
     void selectDT(bool backTrack = false);
-    void changeDT(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 computeSpaceChargeFields(unsigned long long step);
     // void prepareOpalBeamlineSections();
     void dumpStats(long long step, bool psDump, bool statDump);
     void setOptionalVariables();
     bool hasEndOfLineReached(const BoundingBox& globalBoundingBox);
     void handleRestartRun();
 
-    void setTime();
     void doBinaryRepartition();
 
     void transformBunch(const CoordinateSystemTrafo& trafo);
@@ -295,23 +302,70 @@ inline void ParallelTracker::visitTravelingWave(const TravelingWave& as) {
 }
 
 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> x = {Rview(i)[0],Rview(i)[1],Rview(i)[2]};
                              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);
+
+                             // 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();
@@ -332,7 +386,19 @@ inline void ParallelTracker::pushParticles(const BorisPusher& pusher) {
                              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);
-                             pusher.push(x,p,dt);
+
+                                 /** \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
+                                  */
+                             // \TODO check +-
+                             
+                             // pusher.push(x,p,dt);
+                             x = 0.5 * dt * p / Kokkos::sqrt(1.0 + dot(p));
+                             Rview(i) += x;
                          });
 
 
@@ -341,4 +407,4 @@ inline void ParallelTracker::pushParticles(const BorisPusher& pusher) {
     ippl::Comm->barrier();
 }
 
-#endif  // OPAL_ParallelTracker_HH
\ No newline at end of file
+#endif  // OPAL_ParallelTracker_HH