Astra1DDynamic.h

Go to the documentation of this file.

#ifndef OPALX_AstraFIELDMAP1DDYNAMIC_HH
#define OPALX_AstraFIELDMAP1DDYNAMIC_HH
 
#include "Fields/Fieldmap.h"
#include "Physics/Physics.h"
 
#include <Kokkos_Core.hpp>
#include <Kokkos_DualView.hpp>
 
class Astra1DDynamic : public Fieldmap {
public:
    virtual bool getFieldstrength(
            const Vector_t<double, 3>& R, Vector_t<double, 3>& E,
            Vector_t<double, 3>& B) const override;
 
    virtual bool getFieldDerivative(
            const Vector_t<double, 3>& R, Vector_t<double, 3>& E, Vector_t<double, 3>& B,
            const DiffDirection& dir) const override;
 
    virtual void getFieldDimensions(double& zBegin, double& zEnd) const override;
 
    virtual void getFieldDimensions(
            double& xIni, double& xFinal, double& yIni, double& yFinal, double& zIni,
            double& zFinal) const override;
 
    virtual void swap() override;
 
    virtual void getInfo(Inform* msg) override;
 
    virtual double getFrequency() const override;
 
    virtual void setFrequency(double freq) override;
 
    bool isInside(const Vector_t<double, 3>& r) const override {
        return r(2) >= zbegin_m && r(2) < zend_m;
        // && sqrt(r(0) * r(0) + r(1) * r(1)) < rend_m;
    }
 
    template <class ViewType>
    KOKKOS_INLINE_FUNCTION static void computeField(
            const Vector_t<double, 3>& R, Vector_t<double, 3>& E, Vector_t<double, 3>& B,
            const ViewType& FourCoefs, double zbegin, double length, double xlrep, int accuracy) {
        const double RR2 = R(0) * R(0) + R(1) * R(1);
 
        const double two_pi = Physics::two_pi;
        const double pi     = Physics::pi;
        const double dk     = Physics::two_pi / length;
        const double kz     = two_pi * (R(2) - zbegin) / length + pi;
 
        double ez    = FourCoefs(0);
        double ezp   = 0.0;
        double ezpp  = 0.0;
        double ezppp = 0.0;
 
        int n = 1;
        for (int l = 1; l < accuracy; ++l, n += 2) {
            const double base = dk * l;
 
            const double coskzl = Kokkos::cos(kz * l);
            const double sinkzl = Kokkos::sin(kz * l);
 
            ez += FourCoefs(n) * coskzl - FourCoefs(n + 1) * sinkzl;
            ezp += base * (-FourCoefs(n) * sinkzl - FourCoefs(n + 1) * coskzl);
            ezpp += base * base * (-FourCoefs(n) * coskzl + FourCoefs(n + 1) * sinkzl);
            ezppp += base * base * base * (FourCoefs(n) * sinkzl + FourCoefs(n + 1) * coskzl);
        }
 
        const double f  = -(ezpp + ez * xlrep * xlrep) / 16.0;
        const double fp = -(ezppp + ezp * xlrep * xlrep) / 16.0;
 
        const double EfieldR = -(ezp / 2.0 + fp * RR2);
        const double BfieldT = (ez / 2.0 + f * RR2) * xlrep / Physics::c;
 
        E(0) += EfieldR * R(0);
        E(1) += EfieldR * R(1);
        E(2) += ez + 4.0 * f * RR2;
 
        B(0) -= BfieldT * R(1);
        B(1) += BfieldT * R(0);
    }
 
    template <class ViewType>
    KOKKOS_INLINE_FUNCTION static void computeRFField(
            const Vector_t<double, 3>& R, Vector_t<double, 3>& E, Vector_t<double, 3>& B,
            const ViewType& FourCoefs, double zbegin, double zend, double length, double xlrep,
            int accuracy, double electricScale, double magneticScale, double startField,
            double endField) {
        if (R(2) < startField || R(2) >= endField) {
            return;
        }
 
        if (R(2) < zbegin || R(2) >= zend) {
            return;
        }
 
        Vector_t<double, 3> tmpE(0.0), tmpB(0.0);
 
        computeField(R, tmpE, tmpB, FourCoefs, zbegin, length, xlrep, accuracy);
 
        E += electricScale * tmpE;
        B += magneticScale * tmpB;
    }
 
    template <class ViewType>
    KOKKOS_INLINE_FUNCTION static void computeTravelingWaveField(
            const Vector_t<double, 3>& R, Vector_t<double, 3>& E, Vector_t<double, 3>& B,
            const ViewType& FourCoefs, double zbegin, double zend, double length, double xlrep,
            int accuracy, double entryElectricScale, double entryMagneticScale,
            double core1ElectricScale, double core1MagneticScale, double core2ElectricScale,
            double core2MagneticScale, double exitElectricScale, double exitMagneticScale,
            double startCoreField, double startExitField, double mappedStartExitField,
            double periodLength, double cellLength, double elementLength) {
        if (R(2) < -0.5 * periodLength || R(2) + 0.5 * periodLength >= elementLength) {
            return;
        }
 
        Vector_t<double, 3> tmpR({R(0), R(1), R(2) + 0.5 * periodLength});
        Vector_t<double, 3> tmpE(0.0), tmpB(0.0);
 
        if (tmpR(2) < startCoreField) {
            if (!(tmpR(2) >= zbegin && tmpR(2) < zend)) {
                return;
            }
 
            computeField(tmpR, tmpE, tmpB, FourCoefs, zbegin, length, xlrep, accuracy);
            E += entryElectricScale * tmpE;
            B += entryMagneticScale * tmpB;
        } else if (tmpR(2) < startExitField) {
            tmpR(2) -= startCoreField;
            const double z = tmpR(2);
 
            tmpR(2) = tmpR(2) - periodLength * Kokkos::floor(tmpR(2) / periodLength);
            tmpR(2) += startCoreField;
 
            if (!(tmpR(2) >= zbegin && tmpR(2) < zend)) {
                return;
            }
 
            computeField(tmpR, tmpE, tmpB, FourCoefs, zbegin, length, xlrep, accuracy);
            E += core1ElectricScale * tmpE;
            B += core1MagneticScale * tmpB;
 
            tmpE = 0.0;
            tmpB = 0.0;
 
            tmpR(2) = z + cellLength;
            tmpR(2) = tmpR(2) - periodLength * Kokkos::floor(tmpR(2) / periodLength);
            tmpR(2) += startCoreField;
 
            if (!(tmpR(2) >= zbegin && tmpR(2) < zend)) {
                return;
            }
 
            computeField(tmpR, tmpE, tmpB, FourCoefs, zbegin, length, xlrep, accuracy);
            E += core2ElectricScale * tmpE;
            B += core2MagneticScale * tmpB;
        } else {
            tmpR(2) -= mappedStartExitField;
 
            if (!(tmpR(2) >= zbegin && tmpR(2) < zend)) {
                return;
            }
 
            computeField(tmpR, tmpE, tmpB, FourCoefs, zbegin, length, xlrep, accuracy);
            E += exitElectricScale * tmpE;
            B += exitMagneticScale * tmpB;
        }
    }
 
    void applyField(std::shared_ptr<ParticleContainer_t> pc, double) override;
 
    void applyRFField(
            std::shared_ptr<ParticleContainer_t> pc, double electricScale, double magneticScale,
            double startField, double endField);
 
    void applyTravelingWave(
            std::shared_ptr<ParticleContainer_t> pc, double entryElectricScale,
            double entryMagneticScale, double core1ElectricScale, double core1MagneticScale,
            double core2ElectricScale, double core2MagneticScale, double exitElectricScale,
            double exitMagneticScale, double startField, double startCoreField,
            double startExitField, double mappedStartExitField, double periodLength,
            double cellLength, double elementLength);
 
    virtual void getOnaxisEz(std::vector<std::pair<double, double>>& F) override;
 
private:
    Astra1DDynamic(const std::string& filename);
    ~Astra1DDynamic();
 
    void readMap() override;
    void freeMap() override;
 
    Kokkos::DualView<double*> FourCoefs_m;
 
    double frequency_m;
 
    double zbegin_m; 
    double zend_m;   
    double xlrep_m;
 
    double length_m;
 
    int accuracy_m;
 
    int num_gridpz_m;
 
    friend class Fieldmap;
};
 
#endif