537 update multipolet and variablerfcavity to conform with opalx architecture

src/AbsBeamline/ElementBase.cpp

@@ -35,6 +35,7 @@ const std::map<ElementType, std::string> ElementBase::elementTypeToString_s = {
         {ElementType::MARKER, "Marker"},
         {ElementType::MONITOR, "Monitor"},
         {ElementType::MULTIPOLE, "Multipole"},
+        {ElementType::MULTIPOLET, "MULTIPOLET"},
         {ElementType::RFCAVITY, "RFCavity"},
         {ElementType::TRAVELINGWAVE, "TravelingWave"},
         {ElementType::SBEND, "SBEND"},

src/AbsBeamline/MultipoleT.cpp

@@ -20,6 +20,9 @@
 #include "MultipoleTCurvedConstRadius.h"
 #include "MultipoleTStraight.h"
 
+#include "BeamlineGeometry/Geometry.h"
+#include "PartBunch/PartBunch.h"
+
 MultipoleT::MultipoleT(const std::string& name) : ElementBase(name) { chooseImplementation(); }
 
 MultipoleT::MultipoleT(const MultipoleT& right)
@@ -38,6 +41,8 @@ void MultipoleT::accept(BeamlineVisitor& visitor) const {
     visitor.visitMultipoleT(*this);
 }
 
+ElementType MultipoleT::getType() const { return ElementType::MULTIPOLET; }
+
 double MultipoleT::getScaling(const double t) const {
     double scaling = 1.0;
     if (scalingTD_m) {
@@ -159,6 +164,84 @@ Geometry& MultipoleT::getGeometry() { return *implementation_->getGeometry(); }
 
 const Geometry& MultipoleT::getGeometry() const { return *implementation_->getGeometry(); }
 
+Vector_t<double, 3> MultipoleT::bendCoords(const Vector_t<double, 3>& r) const {
+    // The stored frame is the design-orbit entrance tangent, so the arc coordinate is
+    // measured directly. A straight MultipoleT has zero curvature and r comes back unchanged.
+    return GeometryHelper::toBendArcCoords(
+            r, getGeometry().getCurvature(), getGeometry().getElementLength());
+}
+
+bool MultipoleT::isInside(const Vector_t<double, 3>& r) const {
+    const Vector_t<double, 3> arc = bendCoords(r);
+    double zBegin                 = 0.0;
+    double zEnd                   = 0.0;
+    getFieldExtent(zBegin, zEnd);
+    return arc(2) >= zBegin && arc(2) < zEnd && ApertureHelper::isInsideAperture(arc, aperture_m);
+}
+
+size_t MultipoleT::markOutsideAperture(const std::shared_ptr<ParticleContainer_t>& pc) {
+    if (!pc || !getFlagDeleteOnTransverseExit()) {
+        return 0;
+    }
+    const size_t nLocal = pc->getLocalNum();
+    if (nLocal == 0) {
+        return 0;
+    }
+
+    // Members copied to locals; the device kernel must not capture `this`.
+    const ApertureType type = aperture_m.first;
+    const double xLimit     = aperture_m.second[0];
+    const double yLimit     = aperture_m.second[1];
+    const double curvature  = getGeometry().getCurvature();
+    const double bodyLength = getGeometry().getElementLength();
+
+    auto Rview   = pc->R.getView();
+    auto invalid = pc->InvalidMask.getView();
+
+    size_t localMarked = 0;
+    Kokkos::parallel_reduce(
+            "MultipoleT::markOutsideAperture", nLocal,
+            KOKKOS_LAMBDA(const size_t i, size_t& count) {
+                // The aperture belongs to the body, so the window is the geometric body
+                // [0, L) in arc coordinates, not the wider field extent used by isInside().
+                const Vector_t<double, 3> arc =
+                        GeometryHelper::toBendArcCoords(Rview(i), curvature, bodyLength);
+
+                const bool inZ = arc(2) >= 0.0 && arc(2) < bodyLength;
+                const bool hit =
+                        inZ
+                        && !ApertureHelper::isInsideAperture(arc(0), arc(1), type, xLimit, yLimit);
+                const bool newlyMarked = hit && !invalid(i);
+                invalid(i)             = invalid(i) || hit;
+                count += newlyMarked ? 1 : 0;
+            },
+            localMarked);
+    Kokkos::fence();
+
+    return localMarked;
+}
+
+bool MultipoleT::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) {
+    const Vector_t<double, 3> arc = bendCoords(R);
+    double zBegin                 = 0.0;
+    double zEnd                   = 0.0;
+    getFieldExtent(zBegin, zEnd);
+    if (arc(2) < zBegin || arc(2) >= zEnd) {
+        return false;
+    }
+    if (!ApertureHelper::isInsideAperture(arc, aperture_m)) {
+        return true;
+    }
+
+    // The reference particle has to see the same field as every other particle. Without this
+    // the orbit threader walks a straight line through a curved magnet, so the magnet never
+    // stops being the selected element and the threader never reaches the end of the line.
+    apply(R, P, t, E, B);
+    return false;
+}
+
 void MultipoleT::validateConfiguration() const {
     if (2 * config_m.maxFOrder_m + 1 > MultipoleTBase::MaxDerivatives) {
         throw OpalException(

src/AbsBeamline/MultipoleT.h

@@ -1,6 +1,4 @@
 //
-// Cubic Spline Interpolation to replace GSL spline
-//
 // Copyright (c) 2023, Paul Scherrer Institute, Villigen PSI, Switzerland
 // All rights reserved
 //
@@ -82,6 +80,8 @@ class MultipoleT : public ElementBase {
     ElementBase* clone() const override;
     /** Accept a beamline visitor */
     void accept(BeamlineVisitor& visitor) const override;
+    /** @return The element type (ElementType::MULTIPOLET). */
+    ElementType getType() const override;
     /** Return the cell geometry */
     Geometry& getGeometry() override;
     /** Return the cell geometry */
@@ -173,6 +173,10 @@ class MultipoleT : public ElementBase {
     void setElementLength(double length);
     /** Get the length of the magnet */
     double getLength() const { return config_m.length_m; }
+    /** Make the generic ElementBase aperture setter visible alongside the MultipoleT one,
+     *  which would otherwise hide it by name.
+     */
+    using ElementBase::setAperture;
     /** Set the aperture dimensions \n
      * This element only supports a rectangular aperture
      * \param vertAp -> Vertical aperture length
@@ -201,15 +205,32 @@ class MultipoleT : public ElementBase {
     void setBoundingBoxLength(double boundingBoxLength);
     /** Return the longitudinal field-support extent.
      *
-     *  For the current MultipoleT implementation the field support is defined
-     *  on the full local body interval
-     *  latexmath:[z \in [0, L)].
+     *  The body interval latexmath:[[0, L]] plus the reach of the tanh fringe past each
+     *  pole face, measured along the design arc. A hard edge (zero fringe lengths) gives
+     *  the plain body extent.
      */
     void getFieldExtent(double& zBegin, double& zEnd) const override {
-        zBegin = 0.0;
-        zEnd   = getGeometry().getElementLength();
+        zBegin = -MultipoleTBase::FringeReach * config_m.fringeLambdaLeft_m;
+        zEnd   = getGeometry().getElementLength()
+               + MultipoleTBase::FringeReach * config_m.fringeLambdaRight_m;
     }
 
+    /** Return true if the point is inside the field support and the aperture.
+     *
+     *  Both are measured along the design arc, so a curved MultipoleT stays selected as
+     *  the orbit curves through it, the same way SBend does.
+     */
+    bool isInside(const Vector_t<double, 3>& r) const override;
+    /** Mark every particle that leaves the aperture inside the body, in arc coordinates */
+    size_t markOutsideAperture(const std::shared_ptr<ParticleContainer_t>& pc) override;
+    /** Apply the field to the reference particle.
+     *
+     *  @return true if the particle is outside the aperture inside the field extent.
+     */
+    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;
+
     void setScalingName(const std::string& name);
     std::string getScalingName() const { return scalingName_m; }
     void initialiseTimeDependencies() const;
@@ -221,6 +242,8 @@ class MultipoleT : public ElementBase {
     MultipoleTConfig config_m;
 
     void chooseImplementation();
+    /** Map an entrance-frame point to arc coordinates (radial offset, y, arc length) */
+    Vector_t<double, 3> bendCoords(const Vector_t<double, 3>& r) const;
     double getScaling(double t) const;
     void validateConfiguration() const;
 

src/AbsBeamline/MultipoleTBase.h

@@ -101,6 +101,8 @@ class MultipoleTBase {
     static constexpr size_t MaxFactorial         = 20;
     static constexpr size_t MaxPowerInteger      = 20;
     static constexpr unsigned int MaxDerivatives = 20;
+    /// How many fringe lengths past each pole face the tanh fringe is still counted as field.
+    static constexpr double FringeReach = 6.0;
 
     /** Helper function that returns factorial of n for n<=20 on both host and GPU */
     KOKKOS_INLINE_FUNCTION static double factorial(unsigned int n);

src/AbsBeamline/MultipoleTConfig.h

@@ -41,10 +41,14 @@ struct MultipoleTConfig {
     double verticalAperture_m{0.5};
     double horizontalAperture_m{0.5};
     double boundingBoxLength_m{0.0};
-    /** Fringe field parameters **/
-    double fringeS0_m;
-    double fringeLambdaLeft_m;
-    double fringeLambdaRight_m;
+    /** Half separation of the fringe-profile centres [m]; set before field evaluation. */
+    double fringeS0_m{0.0};
+    /** Entrance fringe scale [m]. Zero gives no field-support extension by default;
+     *  the tanh field model requires a positive scale for field evaluation.
+     */
+    double fringeLambdaLeft_m{0.0};
+    /** Exit fringe scale [m], with the same default and requirement as the entrance. */
+    double fringeLambdaRight_m{0.0};
 };
 
-#endif  // OPALX_MULTIPOLETCONFIG_H
+#endif  // OPALX_MULTIPOLETCONFIG_H

src/AbsBeamline/MultipoleTCurvedConstRadius.cpp

@@ -19,11 +19,14 @@
 #include "MultipoleT.h"
 
 MultipoleTCurvedConstRadius::MultipoleTCurvedConstRadius(MultipoleT* element)
-    : MultipoleTBase(element), planarArcGeometry_m(Geometry::makeSBend(1.0, 1.0)) {}
+    : MultipoleTBase(element),
+      planarArcGeometry_m(
+              Geometry::makeSBend(
+                      element->getLength(), element->getBendAngle() / element->getLength())) {}
 
 void MultipoleTCurvedConstRadius::initialise() {
-    planarArcGeometry_m.setElementLength(element_m->getLength());
-    planarArcGeometry_m.setCurvature(element_m->getBendAngle() / element_m->getLength());
+    planarArcGeometry_m = Geometry::makeSBend(
+            element_m->getLength(), element_m->getBendAngle() / element_m->getLength());
     generateTanhCoefficients(element_m->getMaxFOrder() * 2 + 1);
 }
 

src/AbsBeamline/MultipoleTCurvedConstRadius.h

@@ -85,7 +85,7 @@ class MultipoleTCurvedConstRadius final : public MultipoleTBase {
 
 private:
     /** Geometry */
-    Geometry planarArcGeometry_m{Geometry::makeSBend(1.0, 1.0)};
+    Geometry planarArcGeometry_m;
 
     // Helpers
     KOKKOS_INLINE_FUNCTION static Vector_t<double, 3> toMagnetCoords(
@@ -99,15 +99,14 @@ class MultipoleTCurvedConstRadius final : public MultipoleTBase {
 KOKKOS_INLINE_FUNCTION Vector_t<double, 3> MultipoleTCurvedConstRadius::toMagnetCoords(
         const Vector_t<double, 3>& R, const MultipoleTConfig& config) {
     // Skew and entry angle
-    Vector_t<double, 3> result = rotateFrame(R, config);
-    // Go to local Frenet-Serret coordinates
+    const Vector_t<double, 3> rotated = rotateFrame(R, config);
+    // Go to the arc coordinates of the design orbit. This is the same mapping the bends use,
+    // so a positive bend angle turns the orbit towards -x, with the centre of curvature on
+    // the -x side, exactly as for SBEND.
     // Note: if the bend angle is zero, this object is not constructed
-    const double radius  = config.length_m / config.bendAngle_m;
-    const double rMinusX = radius - R[0];
-    const double alpha   = Kokkos::hypot(rMinusX, R[2]);
-    result[0]            = alpha - radius;
-    result[1]            = R[1];
-    result[2]            = radius * Kokkos::atan2(R[2], rMinusX);
+    const double curvature = config.bendAngle_m / config.length_m;
+    Vector_t<double, 3> result =
+            GeometryHelper::toBendArcCoords(rotated, curvature, config.length_m);
     // Magnet origin at the center rather than entry
     result[2] -= config.length_m / 2.0;
     return result;
@@ -161,7 +160,11 @@ KOKKOS_INLINE_FUNCTION bool MultipoleTCurvedConstRadius::computeBField(
             myB[1] += innerSumZ * zzk;
             myB[2] += innerSumS * xszk;
         }
-        B += myB * scaling;
+        // myB is expressed in the basis tangent to the design arc; rotate it into the
+        // element entrance frame, which is the frame the tracker works in.
+        const double sFromEntrance = RPrime[2] + config.length_m / 2.0;
+        B += GeometryHelper::rotateArcFieldToEntry(myB, sFromEntrance, 1.0 / rho, config.length_m)
+             * scaling;
     }
     return !insideAperture;
 }

src/AbsBeamline/VariableRFCavity.h

@@ -50,7 +50,7 @@ class VariableRFCavity : public ElementBase {
     VariableRFCavity& operator=(const VariableRFCavity& rhs);
 
     /** Destructor */
-    virtual ~VariableRFCavity() = default;
+    ~VariableRFCavity() override = default;
 
     /** Return a deep copy */
     ElementBase* clone() const override;

src/BeamlineCore/CMakeLists.txt

@@ -30,6 +30,8 @@ set (HDRS
     TravelingWaveRep.h
     SolenoidRep.h
     ProbeRep.h
+    MultipoleTRep.h
+    VariableRfCavityRep.h
 )
 
 target_include_directories(opalx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})

src/BeamlineCore/MultipoleTRep.h

@@ -0,0 +1,23 @@
+//
+// Copyright (c) 2023, Paul Scherrer Institute, Villigen PSI, Switzerland
+// 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 License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
+
+#ifndef OPALX_MULTIPOLETREP_H
+#define OPALX_MULTIPOLETREP_H
+
+#include "AbsBeamline/MultipoleT.h"
+
+using MultipoleTRep = MultipoleT;
+
+#endif  // OPALX_MULTIPOLETREP_H

src/BeamlineCore/VariableRFCavityRep.h

@@ -0,0 +1,23 @@
+//
+// Copyright (c) 2023, Paul Scherrer Institute, Villigen PSI, Switzerland
+// 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 License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
+
+#ifndef OPALX_VARIABLERFCAVITYREP_H
+#define OPALX_VARIABLERFCAVITYREP_H
+
+#include "AbsBeamline/VariableRFCavity.h"
+
+using VariableRfCavityRep = VariableRFCavity;
+
+#endif  // OPALX_VARIABLERFCAVITYREP_H

src/BeamlineGeometry/Geometry.h

@@ -65,6 +65,10 @@ class Geometry {
     ///@}
 
     GeometryKind kind() const { return kind_m; }
+    /// True for a body that turns the design orbit (sector or rectangular bend).
+    bool isBend() const { return kind_m == GeometryKind::SBend || kind_m == GeometryKind::RBend; }
+    /// True for a rectangular (straight box) bend body.
+    bool isRectangularBend() const { return kind_m == GeometryKind::RBend; }
 
     /// @name Lengths
     ///@{

src/Elements/OpalBeamline.cpp

@@ -224,7 +224,7 @@ void OpalBeamline::save3DLattice() {
 
         mesh.add(*(element.get()));
 
-        if (element->getType() == ElementType::SBEND || element->getType() == ElementType::RBEND) {
+        if (element->getGeometry().isBend()) {
             std::vector<Vector_t<double, 3>> designPath = element->getGeometry().getDesignPath();
             unsigned int size                           = designPath.size();
 

src/Elements/OpalMultipoleT.cpp

@@ -15,8 +15,8 @@
 
 #include "Elements/OpalMultipoleT.h"
 #include <vector>
-#include "AbsBeamline/MultipoleT.h"
 #include "Attributes/Attributes.h"
+#include "BeamlineCore/MultipoleTRep.h"
 
 OpalMultipoleT::OpalMultipoleT()
     : OpalElement(
@@ -103,7 +103,13 @@ void OpalMultipoleT::update() {
     auto tp = Attributes::getRealArray(itsAttr[TP]);
     // Set the attributes
     const auto length = Attributes::getReal(itsAttr[LENGTH]);
-    auto* multT       = dynamic_cast<MultipoleT*>(getElement());
+    if (bendAngle != 0.0 && length == 0.0) {
+        throw OpalException(
+                "OpalMultipoleT::Update",
+                "A curved magnet (non-zero ANGLE) needs a non-zero L: the curvature is "
+                "ANGLE / L, with L the arc length");
+    }
+    auto* multT = dynamic_cast<MultipoleTRep*>(getElement());
     multT->setElementLength(length);
     multT->setBendAngle(bendAngle, varRadius);
     multT->setAperture(Attributes::getReal(itsAttr[VAPERT]), Attributes::getReal(itsAttr[HAPERT]));