Merge branch '349-check-that-element-length-is-long-relative-to-a-time-step-2' into 'master'

Resolve "Check that element length is long relative to a time step"

Closes #349

See merge request !430

src/Algorithms/IndexMap.cpp

@@ -393,9 +393,9 @@ std::pair<double, double> IndexMap::getRange(const IndexMap::value_t::value_type
     return range;
 }
 
-IndexMap::value_t IndexMap::getTouchingElements(const std::pair<double, double> &range) {
-    map_t::iterator it = mapRange2Element_m.begin();
-    const map_t::iterator end = mapRange2Element_m.end();
+IndexMap::value_t IndexMap::getTouchingElements(const std::pair<double, double> &range) const {
+    map_t::const_iterator it = mapRange2Element_m.begin();
+    const map_t::const_iterator end = mapRange2Element_m.end();
     IndexMap::value_t touchingElements;
 
     for (; it != end; ++ it) {

src/Algorithms/IndexMap.h

@@ -54,7 +54,7 @@ public:
     size_t numElements() const;
     std::pair<double, double> getRange(const IndexMap::value_t::value_type &element,
                                        double position) const;
-    IndexMap::value_t getTouchingElements(const std::pair<double, double> &range);
+    IndexMap::value_t getTouchingElements(const std::pair<double, double> &range) const;
 
     class OutOfBounds: public OpalException {
     public:

src/Algorithms/OrbitThreader.cpp

@@ -53,14 +53,15 @@ OrbitThreader::OrbitThreader(const PartData &ref,
                              double maxDiffZBunch,
                              double t,
                              double dt,
-                             double zstop,
+                             StepSizeConfig stepSizes,
                              OpalBeamline &bl) :
     r_m(r),
     p_m(p),
     pathLength_m(s),
     time_m(t),
     dt_m(dt),
-    zstop_m(zstop),
+    stepSizes_m(stepSizes),
+    zstop_m(stepSizes.getFinalZStop() + std::copysign(1.0, dt) * 2 * maxDiffZBunch),
     itsOpalBeamline_m(bl),
     errorFlag_m(0),
     integrator_m(ref),
@@ -105,6 +106,29 @@ OrbitThreader::OrbitThreader(const PartData &ref,
     computeBoundingBox();
 }
 
+void OrbitThreader::checkElementLengths(const std::set<std::shared_ptr<Component>>& fields) {
+    while (!stepSizes_m.reachedEnd() && pathLength_m > stepSizes_m.getZStop()) {
+        ++ stepSizes_m;
+    }
+    if (stepSizes_m.reachedEnd()) {
+        return;
+    }
+    double driftLength = Physics::c * std::abs(stepSizes_m.getdT()) * euclidean_norm(p_m) / Util::getGamma(p_m);
+    for (const std::shared_ptr<Component> field : fields) {
+        double length = field->getElementLength();
+        int numSteps = field->getRequiredNumberOfTimeSteps();
+
+        if (length < numSteps * driftLength) {
+            throw OpalException("OrbitThreader::checkElementLengths",
+                                "The time step is too long compared to the length of the\n"
+                                "element '" + field->getName() + "'\n" +
+                                "The length of the element is: " + std::to_string(length) + "\n"
+                                "The distance the particles drift in " + std::to_string(numSteps) +
+                                " time step(s) is: " + std::to_string(numSteps * driftLength));
+        }
+    }
+}
+
 void OrbitThreader::execute() {
     double initialPathLength = pathLength_m;
 
@@ -123,6 +147,7 @@ void OrbitThreader::execute() {
     auto elementSet = itsOpalBeamline_m.getElements(nextR);
     errorFlag_m = EVERYTHINGFINE;
     do {
+        checkElementLengths(elementSet);
         if (containsCavity(elementSet)) {
             autophaseCavities(elementSet, visitedElements);
         }
@@ -447,4 +472,4 @@ double OrbitThreader::computeDriftLengthToBoundingBox(const std::set<std::shared
     }
 
     return std::numeric_limits<double>::max();
-}
\ No newline at end of file
+}

src/Algorithms/OrbitThreader.h

@@ -23,6 +23,7 @@
 #define OPAL_ORBITTHREADER_H
 
 #include "Algorithms/IndexMap.h"
+#include "Algorithms/StepSizeConfig.h"
 #include "Algorithms/Vektor.h"
 #include "Elements/OpalBeamline.h"
 #include "Steppers/BorisPusher.h"
@@ -41,7 +42,7 @@ public:
                   double maxDiffZBunch,
                   double t,
                   double dT,
-                  double zstop,
+                  StepSizeConfig stepSizes,
                   OpalBeamline &bl);
 
     void execute();
@@ -51,7 +52,7 @@ public:
 
     std::pair<double, double> getRange(const IndexMap::value_t::value_type &element,
                                        double position) const;
-    IndexMap::value_t getTouchingElements(const std::pair<double, double> &range);
+    IndexMap::value_t getTouchingElements(const std::pair<double, double> &range) const;
 
 private:
     /// position of reference particle in lab coordinates
@@ -69,6 +70,7 @@ private:
     double dt_m;
 
     /// final position in path length
+    StepSizeConfig stepSizes_m;
     const double zstop_m;
 
     OpalBeamline &itsOpalBeamline_m;
@@ -112,6 +114,8 @@ private:
     double computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,
                                            const Vector_t & position,
                                            const Vector_t & direction) const;
+
+    void checkElementLengths(const std::set<std::shared_ptr<Component>>& elements);
 };
 
 inline
@@ -127,7 +131,7 @@ std::pair<double, double> OrbitThreader::getRange(const IndexMap::value_t::value
 }
 
 inline
-IndexMap::value_t OrbitThreader::getTouchingElements(const std::pair<double, double> &range) {
+IndexMap::value_t OrbitThreader::getTouchingElements(const std::pair<double, double> &range) const {
     return imap_m.getTouchingElements(range);
 }
 

src/Algorithms/ParallelTTracker.cpp

@@ -247,7 +247,7 @@ void ParallelTTracker::execute() {
                       -rmin(2),
                       itsBunch_m->getT(),
                       (back_track? -minTimeStep: minTimeStep),
-                      stepSizes_m.getFinalZStop() + (back_track? -1: 1) * 2 * rmax(2),
+                      stepSizes_m,
                       itsOpalBeamline_m);
 
     oth.execute();

src/Classic/AbsBeamline/Bend2D.cpp

@@ -208,7 +208,7 @@ void Bend2D::initialise(PartBunchBase<double, 3> *bunch,
         double bendAngleY = 0.0;
         calculateRefTrajectory(bendAngleX, bendAngleY);
         print(msg, bendAngleX, bendAngleY);
-
+        setElementLength(endField_m - startField_m);
         // Pass start and end of field to calling function.
         startField = startField_m;
         endField = endField_m;
@@ -1217,7 +1217,7 @@ bool Bend2D::setupBendGeometry(double &startField, double &endField) {
 
 bool Bend2D::setupDefaultFieldMap() {
 
-    if (length_m <= 0.0) {
+    if (getElementLength() <= 0.0) {
         ERRORMSG("If using \"1DPROFILE1-DEFAULT\" field map you must set the "
                  "chord length of the bend using the L attribute in the OPAL "
                  "input file."

src/Classic/AbsBeamline/BendBase.cpp

@@ -10,7 +10,6 @@ BendBase::BendBase():
 
 BendBase::BendBase(const BendBase &right):
     Component(right),
-    length_m(right.length_m),
     chordLength_m(right.chordLength_m),
     angle_m(right.angle_m),
     entranceAngle_m(right.entranceAngle_m),
@@ -27,7 +26,6 @@ BendBase::BendBase(const BendBase &right):
 
 BendBase::BendBase(const std::string &name):
     Component(name),
-    length_m(0.0),
     chordLength_m(0.0),
     angle_m(0.0),
     entranceAngle_m(0.0),

src/Classic/AbsBeamline/BendBase.h

@@ -17,8 +17,6 @@ public:
     /// Indicates that element bends the beam.
     virtual bool bends() const;
 
-    void setLength(double length);
-    double getLength() const;
     double getChordLength() const;
     virtual void setBendAngle(double angle);
     double getBendAngle() const;
@@ -51,7 +49,6 @@ protected:
     /// Calculate beta*gamma from design energy
     double calcBetaGamma() const;
 
-    double length_m;
     double chordLength_m;
     double angle_m;         ///< Bend angle
     double entranceAngle_m; ///< Angle between incoming reference trajectory
@@ -81,17 +78,6 @@ bool BendBase::bends() const {
     return true;
 }
 
-inline
-void BendBase::setLength(double length) {
-    length_m = std::abs(length);
-}
-
-inline
-double BendBase::getLength() const
-{
-    return length_m;
-}
-
 inline
 double BendBase::getChordLength() const {
     return chordLength_m;

src/Classic/AbsBeamline/Drift.h

@@ -60,6 +60,7 @@ public:
     //set number of slices for map tracking
     std::size_t getNSlices() const; // Philippe was here
 
+    virtual int getRequiredNumberOfTimeSteps() const override;
 private:
 
     double startField_m;
@@ -69,4 +70,10 @@ private:
     void operator=(const Drift &);
 };
 
-#endif // CLASSIC_Drift_HH
+inline
+int Drift::getRequiredNumberOfTimeSteps() const
+{
+    return 1;
+}
+
+#endif // CLASSIC_Drift_HH
\ No newline at end of file

src/Classic/AbsBeamline/ElementBase.h

@@ -372,6 +372,8 @@ public:
 
     virtual BoundingBox getBoundingBoxInLabCoords() const;
 
+    virtual int getRequiredNumberOfTimeSteps() const;
+
 protected:
     bool isInsideTransverse(const Vector_t &r) const;
 
@@ -613,4 +615,10 @@ inline
 bool ElementBase::isElementPositionSet() const
 { return elemedgeSet_m; }
 
+inline
+int ElementBase::getRequiredNumberOfTimeSteps() const
+{
+    return 10;
+}
+
 #endif // CLASSIC_ElementBase_HH
\ No newline at end of file

src/Classic/AbsBeamline/Marker.h

@@ -53,10 +53,18 @@ public:
 
     virtual void getDimensions(double &zBegin, double &zEnd) const override;
 
+    virtual int getRequiredNumberOfTimeSteps() const override;
+
 private:
 
     // Not implemented.
     void operator=(const Marker &);
 };
 
-#endif // CLASSIC_Marker_HH
+inline
+int Marker::getRequiredNumberOfTimeSteps() const
+{
+    return 1;
+}
+
+#endif // CLASSIC_Marker_HH
\ No newline at end of file

src/Classic/AbsBeamline/Monitor.h

@@ -105,6 +105,8 @@ public:
     void setType(Type type);
 
     static void writeStatistics();
+
+    virtual int getRequiredNumberOfTimeSteps() const override;
 private:
 
     // Not implemented.
@@ -125,4 +127,10 @@ void Monitor::setType(Monitor::Type type) {
     type_m = type;
 }
 
+inline
+int Monitor::getRequiredNumberOfTimeSteps() const
+{
+    return 1;
+}
+
 #endif // CLASSIC_Monitor_HH
\ No newline at end of file

src/Classic/AbsBeamline/RBend.cpp

@@ -135,14 +135,14 @@ bool RBend::findChordLength(double &chordLength) {
     const double angle = getBendAngle();
     if (std::abs(angle) > 0.0) {
         double E1 = std::copysign(1.0, angle) * getEntranceAngle();
-        chordLength = 2 * getLength() * std::sin(0.5 * std::abs(angle)) /
+        chordLength = 2 * getElementLength() * std::sin(0.5 * std::abs(angle)) /
             (std::sin(E1) + std::sin(std::abs(angle) - E1));
     } else {
         double refCharge = RefPartBunch_m->getQ();
         double amplitude = (fieldAmplitudeY_m != 0.0) ? fieldAmplitudeY_m : fieldAmplitudeX_m;
         double fieldAmplitude = std::copysign(1.0, refCharge) * std::abs(amplitude);
         double designRadius = calcDesignRadius(fieldAmplitude);
-        chordLength = std::sin(0.5 * (std::asin(getLength() / designRadius - std::sin(getEntranceAngle())) + getEntranceAngle())) * 2 * designRadius;
+        chordLength = std::sin(0.5 * (std::asin(getElementLength() / designRadius - std::sin(getEntranceAngle())) + getEntranceAngle())) * 2 * designRadius;
     }
 
     return true;

src/Classic/AbsBeamline/RBend3D.cpp

@@ -104,7 +104,7 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
         double zBegin = 0.0, zEnd = 0.0;
         fieldmap_m->getFieldDimensions(zBegin, zEnd);
 
-        if (length_m == 0.0) {
+        if (getElementLength() == 0.0) {
             chordLength_m = 0.0;
             double fieldLength = zEnd - zBegin;
             double z = 0.0, dz = fieldLength / 1000;
@@ -121,9 +121,9 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
                 z -= dz;
             }
             chordLength_m = z - zEntryEdge;
-            length_m = chordLength_m;
+            setElementLength(chordLength_m);
         } else {
-            chordLength_m = length_m;
+            chordLength_m = getElementLength();
         }
 
         startField_m = zBegin;
@@ -200,7 +200,7 @@ void RBend3D::goOffline() {
 
 void RBend3D::getDimensions(double &zBegin, double &zEnd) const {
     zBegin = startField_m;
-    zEnd = startField_m + length_m;
+    zEnd = startField_m + getElementLength();
 }
 
 

src/Classic/AbsBeamline/RFCavity.cpp

@@ -61,7 +61,6 @@ RFCavity::RFCavity(const RFCavity &right):
     autophaseVeto_m(right.autophaseVeto_m),
     designEnergy_m(right.designEnergy_m),
     fieldmap_m(right.fieldmap_m),
-    length_m(right.length_m),
     startField_m(right.startField_m),
     endField_m(right.endField_m),
     type_m(right.type_m),
@@ -95,7 +94,6 @@ RFCavity::RFCavity(const std::string &name):
     autophaseVeto_m(false),
     designEnergy_m(-1.0),
     fieldmap_m(nullptr),
-    length_m(0.0),
     startField_m(0.0),
     endField_m(0.0),
     type_m(SW),
@@ -131,7 +129,7 @@ bool RFCavity::apply(const Vector_t &R,
                      Vector_t &E,
                      Vector_t &B) {
     if (R(2) >= startField_m &&
-        R(2) < startField_m + length_m) {
+        R(2) < startField_m + getElementLength()) {
         Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
 
         bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
@@ -151,7 +149,7 @@ bool RFCavity::applyToReferenceParticle(const Vector_t &R,
                                         Vector_t &B) {
 
     if (R(2) >= startField_m &&
-        R(2) < startField_m + length_m) {
+        R(2) < startField_m + getElementLength()) {
         Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
 
         bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
@@ -200,7 +198,7 @@ void RFCavity::initialise(PartBunchBase<double, 3> *bunch, double &startField, d
         }
         frequency_m = fieldmap_m->getFrequency();
     }
-    length_m = endField - startField_m;
+    setElementLength(endField - startField_m);
 }
 
 // In current version ,this function reads in the cavity voltage profile data from file.
@@ -674,7 +672,7 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
     BorisPusher integrator(*RefPartBunch_m->getReference());
     const double cdt = Physics::c * dt;
     const double zbegin = startField_m;
-    const double zend = length_m + startField_m;
+    const double zend = getElementLength() + startField_m;
 
     Vector_t z(0.0, 0.0, zbegin);
     double dz = 0.5 * p(2) / Util::getGamma(p) * cdt;

src/Classic/AbsBeamline/RFCavity.h

@@ -226,7 +226,6 @@ protected:
     double designEnergy_m;
 
     Fieldmap* fieldmap_m;
-    double length_m;
     double startField_m;         /**< starting point of field(m)*/
 
 private:
@@ -510,9 +509,9 @@ CoordinateSystemTrafo RFCavity::getEdgeToBegin() const {
 
 inline
 CoordinateSystemTrafo RFCavity::getEdgeToEnd() const {
-    CoordinateSystemTrafo ret(Vector_t(0, 0, startField_m + length_m),
+    CoordinateSystemTrafo ret(Vector_t(0, 0, startField_m + getElementLength()),
                               Quaternion(1, 0, 0, 0));
     return ret;
 }
 
-#endif // CLASSIC_RFCavity_HH
+#endif // CLASSIC_RFCavity_HH
\ No newline at end of file

src/Classic/AbsBeamline/SBend.cpp

@@ -97,7 +97,6 @@ bool SBend::findChordLength(double &chordLength) {
      * Find bend chord length. If this was not set by the user using the
      * L (length) attribute, infer it from the field map.
      */
-    chordLength = getLength();
+    chordLength = getElementLength();
     return true;
-}
-
+}
\ No newline at end of file

src/Classic/AbsBeamline/Solenoid.cpp

@@ -44,7 +44,6 @@ Solenoid::Solenoid(const Solenoid &right):
     scale_m(right.scale_m),
     scaleError_m(right.scaleError_m),
     startField_m(right.startField_m),
-    length_m(right.length_m),
     fast_m(right.fast_m) {
 }
 
@@ -56,7 +55,6 @@ Solenoid::Solenoid(const std::string &name):
     scale_m(1.0),
     scaleError_m(0.0),
     startField_m(0.0),
-    length_m(0.0),
     fast_m(true) {
 }
 
@@ -89,7 +87,7 @@ bool Solenoid::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B)
 
 bool Solenoid::apply(const Vector_t &R, const Vector_t &/*P*/, const  double &/*t*/, Vector_t &/*E*/, Vector_t &B) {
     if (R(2) >= startField_m
-        && R(2) < startField_m + length_m) {
+        && R(2) < startField_m + getElementLength()) {
         Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
 
         const bool outOfBounds = myFieldmap_m->getFieldstrength(R, tmpE, tmpB);
@@ -104,7 +102,7 @@ bool Solenoid::apply(const Vector_t &R, const Vector_t &/*P*/, const  double &/*
 bool Solenoid::applyToReferenceParticle(const Vector_t &R, const Vector_t &/*P*/, const  double &/*t*/, Vector_t &/*E*/, Vector_t &B) {
 
     if (R(2) >= startField_m
-        && R(2) < startField_m + length_m) {
+        && R(2) < startField_m + getElementLength()) {
         Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
 
         const bool outOfBounds = myFieldmap_m->getFieldstrength(R, tmpE, tmpB);
@@ -131,8 +129,8 @@ void Solenoid::initialise(PartBunchBase<double, 3> *bunch, double &startField, d
         myFieldmap_m->getFieldDimensions(zBegin, zEnd);
 
         startField_m = zBegin;
-        length_m = zEnd - zBegin;
-        endField = startField + length_m;
+        setElementLength(zEnd - zBegin);
+        endField = startField + getElementLength();
     } else {
         endField = startField;
     }
@@ -166,7 +164,7 @@ void Solenoid::setDKS(double ks) {
 
 void Solenoid::getDimensions(double &zBegin, double &zEnd) const {
     zBegin = startField_m;
-    zEnd = startField_m + length_m;
+    zEnd = startField_m + getElementLength();
 }
 
 
@@ -179,13 +177,8 @@ bool Solenoid::isInside(const Vector_t &r) const {
         && myFieldmap_m->isInside(r);
 }
 
-
-double Solenoid::getElementLength() const {
-    return length_m;
-}
-
 void Solenoid::getElementDimensions(double &begin,
                                          double &end) const {
     begin = startField_m;
-    end = begin + length_m;
+    end = begin + getElementLength();
 }
\ No newline at end of file

src/Classic/AbsBeamline/Solenoid.h

@@ -82,8 +82,6 @@ public:
 
     virtual bool isInside(const Vector_t &r) const override;
 
-    virtual double getElementLength() const override;
-
     virtual void getElementDimensions(double &zBegin, double &zEnd) const override;
 
     virtual CoordinateSystemTrafo getEdgeToBegin() const override;
@@ -98,7 +96,6 @@ private:
     double scaleError_m;                /**< scale multiplier error*/
 
     double startField_m;           /**< startingpoint of field, m*/
-    double length_m;
 
     bool fast_m;
     // Not implemented.
@@ -119,9 +116,9 @@ CoordinateSystemTrafo Solenoid::getEdgeToBegin() const
 inline
 CoordinateSystemTrafo Solenoid::getEdgeToEnd() const
 {
-    CoordinateSystemTrafo ret(Vector_t(0, 0, startField_m + length_m),
+    CoordinateSystemTrafo ret(Vector_t(0, 0, startField_m + getElementLength()),
                               Quaternion(1, 0, 0, 0));
 
     return ret;
 }
-#endif // CLASSIC_Solenoid_HH
+#endif // CLASSIC_Solenoid_HH
\ No newline at end of file

src/Classic/AbsBeamline/Source.h

@@ -39,6 +39,8 @@ public:
 
     virtual void getDimensions(double &zBegin, double &zEnd) const override;
 
+    virtual int getRequiredNumberOfTimeSteps() const override;
+
     void setTransparent();
 private:
 
@@ -52,4 +54,10 @@ private:
     // Not implemented.
     void operator=(const Source &);
 };
+
+inline
+int Source::getRequiredNumberOfTimeSteps() const
+{
+    return 0;
+}
 #endif // CLASSIC_SOURCE_HH
\ No newline at end of file