field.hpp

Go to the documentation of this file.

#ifndef FIELD_HPP
#define FIELD_HPP
#include "space_settings.hpp"
#include "NamelistReader.hpp"

#include "globals.hpp"
#include "gyro_radius.hpp"

/* VarType specifies whether the field is a vector, a 2D vector, or a scalar
 * */
enum class VarType{
    Vector,
    Vector2D,
    Scalar
};

/* Compile-time function to use 2D vector if field-following planes are not used (i.e. XGCa) */
template<PhiInterpType PIT> constexpr VarType vec2d_if_axisym();
template<> constexpr VarType vec2d_if_axisym<PhiInterpType::Planes>(){return VarType::Vector;}
template<> constexpr VarType vec2d_if_axisym<PhiInterpType::None>(){return VarType::Vector2D;}

/* Coordinate correction for field interpolations */
class FieldCorrection{
    const double gamma_psi;

    public:
    double r0;
    double r1;
    double z0;
    double z1;

    KOKKOS_INLINE_FUNCTION FieldCorrection(const SimdVector2D& gradpsi, int i_simd)
      : gamma_psi(1.0/gradpsi.magnitude(i_simd)) {}

    KOKKOS_INLINE_FUNCTION void set(int i_simd, double basis, const SimdVector2D &gradpsi){
        r0 = basis + (1.0-basis) * gamma_psi *  gradpsi.r[i_simd];
        r1 =         (1.0-basis) * gamma_psi *  gradpsi.z[i_simd];
        z0 =         (1.0-basis) * gamma_psi *(-gradpsi.z[i_simd]);
        z1 = basis + (1.0-basis) * gamma_psi *  gradpsi.r[i_simd];
    }
};


template<VarType VT, PhiInterpType PIT>
struct Field{};

//Struct used for E_phi_ff, dAs_phi_ff, dAh_phi_ff, As_phi_ff, Ah_phi_ff. 3 vector components, 2 contributing planes 
template<>
struct Field<VarType::Vector,PhiInterpType::Planes>{
    double V[3][2];

    static constexpr PhiInterpType PIT = PhiInterpType::Planes;
    static constexpr VarType VT = VarType::Vector;

    static constexpr int NDIM(){return 3;}
    static constexpr int NPHI(){return 2;}
    static constexpr int SIZE(){return NDIM()*NPHI();}

    KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts, const FieldCorrection& corr) const{
        const double& phi_w0 = phi_wts.phi.w[0][i_simd];
        const double& phi_w1 = phi_wts.phi.w[1][i_simd];
        vec.r[i_simd]   += wp*(phi_w0*( corr.r0*V[PIR][0] + corr.z0*V[PIZ][0] )
                             + phi_w1*( corr.r0*V[PIR][1] + corr.z0*V[PIZ][1] ) );
        vec.z[i_simd]   += wp*(phi_w0*( corr.r1*V[PIR][0] + corr.z1*V[PIZ][0] )
                             + phi_w1*( corr.r1*V[PIR][1] + corr.z1*V[PIZ][1] ) );
        vec.phi[i_simd] += wp*(phi_w0*  V[PIP][0]
                             + phi_w1*  V[PIP][1] );
    }

    KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const SimdGyroWeights<GyroKin>& rho_wts, const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts, const FieldCorrection& corr, const Field<VarType::Vector,PhiInterpType::Planes>& field2) const{
        const double& phi_w0 = phi_wts.phi.w[0][i_simd];
        const double& phi_w1 = phi_wts.phi.w[1][i_simd];
#ifdef NEWGYROMATRIX
        const double& rho_w00 = rho_wts.rho_wts[0].w[0][i_simd];
        const double& rho_w01 = rho_wts.rho_wts[0].w[1][i_simd];
        const double& rho_w10 = rho_wts.rho_wts[1].w[0][i_simd];
        const double& rho_w11 = rho_wts.rho_wts[1].w[1][i_simd];
#else
        const double& rho_w00 = rho_wts.rho_wts.w[0][i_simd];
        const double& rho_w01 = rho_wts.rho_wts.w[1][i_simd];
        const double& rho_w10 = rho_wts.rho_wts.w[0][i_simd]; // = w00
        const double& rho_w11 = rho_wts.rho_wts.w[1][i_simd]; // = w01
#endif

        vec.r[i_simd]   += wp*(phi_w0*( corr.r0*(rho_w00*V[PIR][0] + rho_w01*field2.V[PIR][0])
                                      + corr.z0*(rho_w00*V[PIZ][0] + rho_w01*field2.V[PIZ][0]) )
                             + phi_w1*( corr.r0*(rho_w10*V[PIR][1] + rho_w11*field2.V[PIR][1])
                                      + corr.z0*(rho_w10*V[PIZ][1] + rho_w11*field2.V[PIZ][1]) ) );
        vec.z[i_simd]   += wp*(phi_w0*( corr.r1*(rho_w00*V[PIR][0] + rho_w01*field2.V[PIR][0])
                                      + corr.z1*(rho_w00*V[PIZ][0] + rho_w01*field2.V[PIZ][0]) )
                             + phi_w1*( corr.r1*(rho_w10*V[PIR][1] + rho_w11*field2.V[PIR][1])
                                      + corr.z1*(rho_w10*V[PIZ][1] + rho_w11*field2.V[PIZ][1]) ) );
        vec.phi[i_simd] += wp*(phi_w0*  (rho_w00*V[PIP][0] + rho_w01*field2.V[PIP][0])
                             + phi_w1*  (rho_w10*V[PIP][1] + rho_w11*field2.V[PIP][1]) );
    }

};

//Struct used for E00_ff. 2 vector components, 2 contributing planes
template<>
struct Field<VarType::Vector2D,PhiInterpType::Planes>{
    double V[2][2];

    static constexpr PhiInterpType PIT = PhiInterpType::Planes;
    static constexpr VarType VT = VarType::Vector2D;

    static constexpr int NDIM(){return 2;}
    static constexpr int NPHI(){return 2;}
    static constexpr int SIZE(){return NDIM()*NPHI();}

    KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts, const FieldCorrection& corr) const{
        const double& phi_w0 = phi_wts.phi.w[0][i_simd];
        const double& phi_w1 = phi_wts.phi.w[1][i_simd];
        vec.r[i_simd]   += wp*(phi_w0*( corr.r0*V[PIR][0] + corr.z0*V[PIZ][0] )
                             + phi_w1*( corr.r0*V[PIR][1] + corr.z0*V[PIZ][1] ) );
        vec.z[i_simd]   += wp*(phi_w0*( corr.r1*V[PIR][0] + corr.z1*V[PIZ][0] )
                             + phi_w1*( corr.r1*V[PIR][1] + corr.z1*V[PIZ][1] ) );
    }
};

//Struct used for ddpotdt_phi_ff, As_phi_ff, Ah_phi_ff for EM.  2 contributing planes
template<>
struct Field<VarType::Scalar,PhiInterpType::Planes>{
    double S[2];

    static constexpr PhiInterpType PIT = PhiInterpType::Planes;
    static constexpr VarType VT = VarType::Scalar;

    static constexpr int NDIM(){return 1;}
    static constexpr int NPHI(){return 2;}
    static constexpr int SIZE(){return NDIM()*NPHI();}

    KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp, const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts) const{
        const double& phi_w0 = phi_wts.phi.w[0][i_simd];
        const double& phi_w1 = phi_wts.phi.w[1][i_simd];
        sca[i_simd] += wp* ( phi_w0*S[0] + phi_w1*S[1] );
    }

    KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp, const SimdGyroWeights<GyroKin>& rho_wts, const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts, const Field<VarType::Scalar,PhiInterpType::Planes>& field2) const{
        const double& phi_w0 = phi_wts.phi.w[0][i_simd];
        const double& phi_w1 = phi_wts.phi.w[1][i_simd];
#ifdef NEWGYROMATRIX
        const double& rho_w00 = rho_wts.rho_wts[0].w[0][i_simd];
        const double& rho_w01 = rho_wts.rho_wts[0].w[1][i_simd];
        const double& rho_w10 = rho_wts.rho_wts[1].w[0][i_simd];
        const double& rho_w11 = rho_wts.rho_wts[1].w[1][i_simd];
#else
        const double& rho_w00 = rho_wts.rho_wts.w[0][i_simd];
        const double& rho_w01 = rho_wts.rho_wts.w[1][i_simd];
        const double& rho_w10 = rho_wts.rho_wts.w[0][i_simd]; // = w00
        const double& rho_w11 = rho_wts.rho_wts.w[1][i_simd]; // = w01
#endif

        sca[i_simd] += wp* (  phi_w0*rho_w00*S[0]
                            + phi_w1*rho_w10*S[1]
                            + phi_w0*rho_w01*field2.S[0]
                            + phi_w1*rho_w11*field2.S[1]);
    }

    // Scatter
    KOKKOS_INLINE_FUNCTION void scatter(const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts, int i_simd, double wp){
        const double& phi_w0 = phi_wts.phi.w[0][i_simd];
        const double& phi_w1 = phi_wts.phi.w[1][i_simd];
        access_add(&(S[0]), wp*phi_w0);
        access_add(&(S[1]), wp*phi_w1);
    }

    // Scatter (gyrokinetic)
    KOKKOS_INLINE_FUNCTION void scatter(const SimdGyroWeights<GyroKin>& rho_wts, const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts, int i_simd, double wp, Field<VarType::Scalar,PhiInterpType::Planes>& field2){
        const double& phi_w0 = phi_wts.phi.w[0][i_simd];
        const double& phi_w1 = phi_wts.phi.w[1][i_simd];
#ifdef NEWGYROMATRIX
        const double& rho_w00 = rho_wts.rho_wts[0].w[0][i_simd];
        const double& rho_w01 = rho_wts.rho_wts[0].w[1][i_simd];
        const double& rho_w10 = rho_wts.rho_wts[1].w[0][i_simd];
        const double& rho_w11 = rho_wts.rho_wts[1].w[1][i_simd];
#else
        const double& rho_w00 = rho_wts.rho_wts.w[0][i_simd];
        const double& rho_w01 = rho_wts.rho_wts.w[1][i_simd];
        const double& rho_w10 = rho_wts.rho_wts.w[0][i_simd]; // = w00
        const double& rho_w11 = rho_wts.rho_wts.w[1][i_simd]; // = w01
#endif

        access_add(&(S[0]), wp*phi_w0*rho_w00);
        access_add(&(S[1]), wp*phi_w1*rho_w10);
        access_add(&(field2.S[0]), wp*phi_w0*rho_w01);
        access_add(&(field2.S[1]), wp*phi_w1*rho_w11);
    }

    Field<VarType::Scalar,PhiInterpType::Planes>& operator+=(const Field<VarType::Scalar,PhiInterpType::Planes>& rhs){
        S[0] += rhs.S[0];
        S[1] += rhs.S[1];
        return *this;
    }
};

// Struct used for E_rho, dEr_B2, dEz_B2, and du2_E_rho. 2 vector components (r,z)
template<>
struct Field<VarType::Vector2D,PhiInterpType::None>{
    double E[2];

    static constexpr PhiInterpType PIT = PhiInterpType::None;
    static constexpr VarType VT = VarType::Vector2D;

    static constexpr int NDIM(){return 2;}
    static constexpr int NPHI(){return 1;}
    static constexpr int SIZE(){return NDIM()*NPHI();}

    KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const FieldCorrection& corr) const{
        vec.r[i_simd] += wp*( corr.r0*E[PIR] + corr.z0*E[PIZ] );
        vec.z[i_simd] += wp*( corr.r1*E[PIR] + corr.z1*E[PIZ] );
    }

    KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const SimdGyroWeights<GyroKin>& rho_wts, const FieldCorrection& corr, const Field<VarType::Vector2D,PhiInterpType::None>& field2) const{
        const double& rho_w00 = rho_wts.rho_wts.w[0][i_simd];
        const double& rho_w01 = rho_wts.rho_wts.w[1][i_simd];
        vec.r[i_simd] += wp*( corr.r0*(rho_w00*E[PIR] + rho_w01*field2.E[PIR])
                             +corr.z0*(rho_w00*E[PIZ] + rho_w01*field2.E[PIZ]) );
        vec.z[i_simd] += wp*( corr.r1*(rho_w00*E[PIR] + rho_w01*field2.E[PIR])
                             +corr.z1*(rho_w00*E[PIZ] + rho_w01*field2.E[PIZ]) );
    }

};

// Struct for vector without planar interpolation
template<>
struct Field<VarType::Vector,PhiInterpType::None>{
    double V[3];

    static constexpr PhiInterpType PIT = PhiInterpType::None;
    static constexpr VarType VT = VarType::Vector;

    static constexpr int NDIM(){return 3;}
    static constexpr int NPHI(){return 1;}
    static constexpr int SIZE(){return NDIM()*NPHI();}

    KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp) const{
        vec.r[i_simd]   += wp*V[PIR];
        vec.z[i_simd]   += wp*V[PIZ];
        vec.phi[i_simd] += wp*V[PIP];
    }
};

// Struct for scalar without planar interpolation (i.e. just a double)
template<>
struct Field<VarType::Scalar,PhiInterpType::None>{
    double S;

    static constexpr PhiInterpType PIT = PhiInterpType::None;
    static constexpr VarType VT = VarType::Scalar;

    static constexpr int NDIM(){return 1;}
    static constexpr int NPHI(){return 1;}
    static constexpr int SIZE(){return NDIM()*NPHI();}

    KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp){
        sca[i_simd] += wp*S;
    }

    // Scatter
    KOKKOS_INLINE_FUNCTION void scatter(const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts, int i_simd, double wp){
        access_add(&S, wp);
    }

    // Scatter (gyrokinetic)
    KOKKOS_INLINE_FUNCTION void scatter(const SimdGyroWeights<GyroKin>& rho_wts, const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts, int i_simd, double wp, Field<VarType::Scalar,PhiInterpType::None>& field2){
        const double& rho_w00 = rho_wts.rho_wts.w[0][i_simd];
        const double& rho_w01 = rho_wts.rho_wts.w[1][i_simd];
        access_add(&(       S),wp*rho_w00);
        access_add(&(field2.S),wp*rho_w01);
    }

    /* Cleaner access patterns if interacting with the same type */

    KOKKOS_INLINE_FUNCTION void operator += (const Field<VarType::Scalar,PhiInterpType::None>& rhs){
        S += rhs.S;
    }

    /* Cleaner access patterns if interacting with a double */

    // Conversion to double
    KOKKOS_INLINE_FUNCTION explicit operator double() const {
        return S;
    }

    // Assignment from double
    KOKKOS_INLINE_FUNCTION void operator = (const double val) {
        S = val;
    }

    // Add double
    KOKKOS_INLINE_FUNCTION void operator += (const double val){
        S += val;
    }
};

#endif