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]) );
    }
};
 
// Returns false if any value is not finite
KOKKOS_INLINE_FUNCTION bool isfinite(const Field<VarType::Vector,PhiInterpType::Planes>& f){
    return (isfinite(f.V[0][0]) && isfinite(f.V[0][1]) && isfinite(f.V[1][0]) &&
            isfinite(f.V[1][1]) && isfinite(f.V[2][0]) && isfinite(f.V[2][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;
    }
 
    // Atomic += wrapper
/*    KOKKOS_INLINE_FUNCTION void access_add_field(const Field<VarType::Scalar,PhiInterpType::Planes>& rhs){
        access_add(&(S[0]), rhs.S[0]);
        access_add(&(S[1]), rhs.S[1]);
    }
 
    // Default constructor
    KOKKOS_INLINE_FUNCTION Field<VarType::Scalar,PhiInterpType::Planes>(){}
 
    // Constructor set to scalar
    KOKKOS_INLINE_FUNCTION Field<VarType::Scalar,PhiInterpType::Planes>(double initval){
        S[0] = initval;
        S[1] = initval;
    }*/
};
 
// 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) const{
        sca[i_simd] += wp*S;
    }
 
    KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp, const SimdPhiWeights<get_phi_wt_usage(PIT)>& phi_wts) const{
        gather(sca, i_simd, wp);
    }
 
    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,PIT>& field2) const{
        const double& rho_w0 = rho_wts.rho_wts.w[0][i_simd];
        const double& rho_w1 = rho_wts.rho_wts.w[1][i_simd];
 
        sca[i_simd] += wp*(rho_w0*S + rho_w1*field2.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;
    }
 
/*    // Atomic += wrapper
    KOKKOS_INLINE_FUNCTION void access_add_field(const Field<VarType::Scalar,PhiInterpType::None>& rhs){
        access_add(&(S), rhs.S);
    }
 
    // Default constructor
    KOKKOS_INLINE_FUNCTION Field<VarType::Scalar,PhiInterpType::None>(){}
 
    // Constructor set to scalar
    KOKKOS_INLINE_FUNCTION Field<VarType::Scalar,PhiInterpType::None>(double initval){
        S = initval;
    }*/
};
 
typedef Field<VarType::Scalar,PhiInterpType::None> ScalarField;
 
#endif