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 "grid.hpp"
#include "linear_weights.hpp"
#include "local_fields.hpp"
#include "push_controls.hpp"
#include "field_weights.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 double (&wphi)[2], const FieldCorrection& corr) const{
        vec.r[i_simd]   += wp*(wphi[0]*( corr.r0*V[PIR][0] + corr.z0*V[PIZ][0] )
                             + wphi[1]*( corr.r0*V[PIR][1] + corr.z0*V[PIZ][1] ) );
        vec.z[i_simd]   += wp*(wphi[0]*( corr.r1*V[PIR][0] + corr.z1*V[PIZ][0] )
                             + wphi[1]*( corr.r1*V[PIR][1] + corr.z1*V[PIZ][1] ) );
        vec.phi[i_simd] += wp*(wphi[0]*  V[PIP][0]
                             + wphi[1]*  V[PIP][1] );
    }

    KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const FieldWeights<GyroKin, PhiInterpType::Planes>& wts, const FieldCorrection& corr, const Field<VarType::Vector,PhiInterpType::Planes>& field2) const{
# ifdef NEWGYROMATRIX
        vec.r[i_simd]   += wp*(wts.phi.w[0]*( corr.r0*(wts.rho[0].w[0]*V[PIR][0] + wts.rho[0].w[1]*field2.V[PIR][0])
                                            + corr.z0*(wts.rho[0].w[0]*V[PIZ][0] + wts.rho[0].w[1]*field2.V[PIZ][0]) )
                             + wts.phi.w[1]*( corr.r0*(wts.rho[1].w[0]*V[PIR][1] + wts.rho[1].w[1]*field2.V[PIR][1])
                                            + corr.z0*(wts.rho[1].w[0]*V[PIZ][1] + wts.rho[1].w[1]*field2.V[PIZ][1]) ) );
        vec.z[i_simd]   += wp*(wts.phi.w[0]*( corr.r1*(wts.rho[0].w[0]*V[PIR][0] + wts.rho[0].w[1]*field2.V[PIR][0])
                                            + corr.z1*(wts.rho[0].w[0]*V[PIZ][0] + wts.rho[0].w[1]*field2.V[PIZ][0]) )
                             + wts.phi.w[1]*( corr.r1*(wts.rho[1].w[0]*V[PIR][1] + wts.rho[1].w[1]*field2.V[PIR][1])
                                            + corr.z1*(wts.rho[1].w[0]*V[PIZ][1] + wts.rho[1].w[1]*field2.V[PIZ][1]) ) );
        vec.phi[i_simd] += wp*(wts.phi.w[0]*  (wts.rho[0].w[0]*V[PIP][0] + wts.rho[0].w[1]*field2.V[PIP][0])
                             + wts.phi.w[1]*  (wts.rho[1].w[0]*V[PIP][1] + wts.rho[1].w[1]*field2.V[PIP][1]) );
# else
        vec.r[i_simd]   += wp*(wts.phi.w[0]*( corr.r0*(wts.rho.w[0]*V[PIR][0] + wts.rho.w[1]*field2.V[PIR][0])
                                            + corr.z0*(wts.rho.w[0]*V[PIZ][0] + wts.rho.w[1]*field2.V[PIZ][0]) )
                             + wts.phi.w[1]*( corr.r0*(wts.rho.w[0]*V[PIR][1] + wts.rho.w[1]*field2.V[PIR][1])
                                            + corr.z0*(wts.rho.w[0]*V[PIZ][1] + wts.rho.w[1]*field2.V[PIZ][1]) ) );
        vec.z[i_simd]   += wp*(wts.phi.w[0]*( corr.r1*(wts.rho.w[0]*V[PIR][0] + wts.rho.w[1]*field2.V[PIR][0])
                                            + corr.z1*(wts.rho.w[0]*V[PIZ][0] + wts.rho.w[1]*field2.V[PIZ][0]) )
                             + wts.phi.w[1]*( corr.r1*(wts.rho.w[0]*V[PIR][1] + wts.rho.w[1]*field2.V[PIR][1])
                                            + corr.z1*(wts.rho.w[0]*V[PIZ][1] + wts.rho.w[1]*field2.V[PIZ][1]) ) );
        vec.phi[i_simd] += wp*(wts.phi.w[0]*  (wts.rho.w[0]*V[PIP][0] + wts.rho.w[1]*field2.V[PIP][0])
                             + wts.phi.w[1]*  (wts.rho.w[0]*V[PIP][1] + wts.rho.w[1]*field2.V[PIP][1]) );
# endif
    }

};

//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 double (&wphi)[2], const FieldCorrection& corr) const{
        vec.r[i_simd]   += wp*(wphi[0]*( corr.r0*V[PIR][0] + corr.z0*V[PIZ][0] )
                             + wphi[1]*( corr.r0*V[PIR][1] + corr.z0*V[PIZ][1] ) );
        vec.z[i_simd]   += wp*(wphi[0]*( corr.r1*V[PIR][0] + corr.z1*V[PIZ][0] )
                             + wphi[1]*( 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 double (&wphi)[2]) const{
        sca[i_simd] += wp* ( wphi[0]*S[0] + wphi[1]*S[1] );
    }

    KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp, const FieldWeights<GyroKin, PhiInterpType::Planes>& wts, const Field<VarType::Scalar,PhiInterpType::Planes>& field2) const{
# ifdef NEWGYROMATRIX
        sca[i_simd] += wp* (  wts.phi.w[0]*wts.rho[0].w[0]*S[0]
                            + wts.phi.w[1]*wts.rho[1].w[0]*S[1]
                            + wts.phi.w[0]*wts.rho[0].w[1]*field2.S[0]
                            + wts.phi.w[1]*wts.rho[1].w[1]*field2.S[1]);
# else
        sca[i_simd] += wp* (  wts.phi.w[0]*wts.rho.w[0]*S[0]
                            + wts.phi.w[1]*wts.rho.w[0]*S[1]
                            + wts.phi.w[0]*wts.rho.w[1]*field2.S[0]
                            + wts.phi.w[1]*wts.rho.w[1]*field2.S[1]);
# endif
    }

    // Scatter
    KOKKOS_INLINE_FUNCTION void scatter(const FieldWeights<DriftKin, PhiInterpType::Planes>& wts, double wp){
        access_add(&(S[0]), wp*wts.phi.w[0]);
        access_add(&(S[1]), wp*wts.phi.w[1]);
    }

    // Scatter (gyrokinetic)
    KOKKOS_INLINE_FUNCTION void scatter(const FieldWeights<GyroKin, PhiInterpType::Planes>& wts, double wp, Field<VarType::Scalar,PhiInterpType::Planes>& field2){
#ifdef NEWGYROMATRIX
        access_add(&(S[0]), wp*wts.phi.w[0]*wts.rho[0].w[0]);
        access_add(&(S[1]), wp*wts.phi.w[1]*wts.rho[1].w[0]);
        access_add(&(field2.S[0]), wp*wts.phi.w[0]*wts.rho[0].w[1]);
        access_add(&(field2.S[1]), wp*wts.phi.w[1]*wts.rho[1].w[1]);
#else
        access_add(&(S[0]), wp*wts.phi.w[0]*wts.rho.w[0]);
        access_add(&(S[1]), wp*wts.phi.w[1]*wts.rho.w[0]);
        access_add(&(field2.S[0]), wp*wts.phi.w[0]*wts.rho.w[1]);
        access_add(&(field2.S[1]), wp*wts.phi.w[1]*wts.rho.w[1]);
#endif
    }

    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 FieldWeights<GyroKin, PhiInterpType::None>& wts, const FieldCorrection& corr, const Field<VarType::Vector2D,PhiInterpType::None>& field2) const{
        vec.r[i_simd] += wp*( corr.r0*(wts.rho.w[0]*E[PIR] + wts.rho.w[1]*field2.E[PIR])
                             +corr.z0*(wts.rho.w[0]*E[PIZ] + wts.rho.w[1]*field2.E[PIZ]) );
        vec.z[i_simd] += wp*( corr.r1*(wts.rho.w[0]*E[PIR] + wts.rho.w[1]*field2.E[PIR])
                             +corr.z1*(wts.rho.w[0]*E[PIZ] + wts.rho.w[1]*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 FieldWeights<DriftKin, PhiInterpType::None>& wts, double wp){
        access_add(&S, wp);
    }

    // Scatter (gyrokinetic)
    KOKKOS_INLINE_FUNCTION void scatter(const FieldWeights<GyroKin, PhiInterpType::None>& wts, double wp, Field<VarType::Scalar,PhiInterpType::None>& field2){
        access_add(&(S),wp*wts.rho.w[0]);
        access_add(&(field2.S),wp*wts.rho.w[1]);
    }

    /* 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