magnetic_field.hpp

Go to the documentation of this file.

#ifndef MAGNETIC_FIELD_HPP
#define MAGNETIC_FIELD_HPP
#include "magnetic_equil_files.hpp"
#include "equil.hpp"
#include "bicub.hpp"
#include "cub_interp.hpp"
#include "range_view.hpp"
#include "tricub.hpp"

// Magnetic field class
template<class Device>
class MagneticField {
    // This line enables private access between the host and device implementations of this class
    // i.e. Makes the copy from host to device more straightforward
    template<class Device2> friend class MagneticField;

    KOKKOS_INLINE_FUNCTION void derivs(const double (&x)[2],double phi,double (&dx)[2]) const;

#ifndef STELLARATOR
    KOKKOS_INLINE_FUNCTION double I_value(double psi_in,int rgn3) const;
    KOKKOS_INLINE_FUNCTION double I_deriv(double psi_in,int rgn3) const;
#endif

    int ff_step;     
    int ff_order;    

#ifndef STELLARATOR
    Bicub<Device> psi_bicub;   
    CubInterp<Device> I_interp;   
#else
    Tricub<Device> psi_tricub;
    Tricub<Device> Br_tricub;
    Tricub<Device> Bz_tricub;
    Tricub<Device> Bphi_tricub;
#endif

    public:

    double bt_sign;  
    double bp_sign;  
    RZBounds bounds; 
    double inpsi;    
    double outpsi;   
    Equilibrium equil; 

    MagneticField(NLReader::NamelistReader& nlr, const DomainDecomposition<DeviceType>& pol_decomp,
                         const MagneticEquilFiles::Ptr& equil_files, const View<Equil::XPoint*, HostType>& xpts,
                         const RZPair& axis_in, const RZPair& psi_ref_point, const View<double*, HostType>& psiN);

    void write() const;

    // Create a mirror with a different device type
    template<class Device2>
    inline MagneticField<Device2> mirror() const{
        MagneticField<Device2> m;
        m.bt_sign = bt_sign;
        m.bp_sign = bp_sign;
        m.ff_step = ff_step;
        m.ff_order = ff_order;
        m.bounds = bounds;

        m.inpsi = inpsi;
        m.outpsi = outpsi;

        m.equil = equil;
#ifndef STELLARATOR
        m.psi_bicub = psi_bicub.template mirror<Device2>();
        m.I_interp = I_interp.template mirror<Device2>();
#else
        m.psi_tricub = psi_tricub.template mirror<Device2>();
        m.Br_tricub = Br_tricub.template mirror<Device2>();
        m.Bz_tricub = Bz_tricub.template mirror<Device2>();
        m.Bphi_tricub = Bphi_tricub.template mirror<Device2>();
#endif

        return m;
    }

    // Analytic constructors
#ifdef STELLARATOR
    MagneticField(PsiOption psi_opt, int eq_mr_in=-1, int eq_mz_in=-1, int eq_mphi_in=-1);
#else
    MagneticField(PsiOption psi_opt, double safety_factor_coeff=1.0, int eq_mr_in=-1, int eq_mz_in=-1, int eq_mpsi_in=-1);
#endif

    // Default constructor
    MagneticField(){}

    KOKKOS_INLINE_FUNCTION double psi_norm() const;

    KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector& v,Simd<double>& psi_out) const;

    KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector2D& x,const Simd<double>& phi, Simd<double>& psi_out) const;

    KOKKOS_INLINE_FUNCTION double get_psi(double r, double z, double phi) const;

    KOKKOS_INLINE_FUNCTION void get_psi_and_derivs(double r, double z, double phi, double& psi, double& dpsidr, double& dpsidz, double& dpsidphi) const;

    KOKKOS_INLINE_FUNCTION void get_psi_unit_vec(const SimdVector2D& x, double phi, Simd<double>& cosa, Simd<double>& sina) const;

    KOKKOS_INLINE_FUNCTION double geometry_r(double r) const;

    KOKKOS_INLINE_FUNCTION bool is_below_xpt_tangent(double r,double z) const;
    KOKKOS_INLINE_FUNCTION bool is_above_xpt2_tangent(double r, double z) const;
    KOKKOS_INLINE_FUNCTION bool is_in_region_3a_or_sep_leg(double r,double z,double psi) const;
    KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r,double z,double psi) const;
    KOKKOS_INLINE_FUNCTION bool is_in_region_1(double r,double z,double psi) const;
    KOKKOS_INLINE_FUNCTION void check_boundaries(const SimdVector2D& x, Simd<bool>& rz_outside) const;
    KOKKOS_INLINE_FUNCTION void get_theta(const SimdVector2D& x, Simd<double>& theta) const;
    KOKKOS_INLINE_FUNCTION RZPair get_xpoint(int ixpt) const;
    KOKKOS_INLINE_FUNCTION RZPair get_axis() const;

    // Get complete B-field info
    KOKKOS_INLINE_FUNCTION void field(const SimdVector& v, SimdVector &bvec, SimdVector (&jacb)[3], Simd<double>& psivec, SimdVector2D &gradpsi, Simd<bool>& rz_outside) const;

    KOKKOS_INLINE_FUNCTION void bvec_interpol(double r,double z,double phi,double &br,double &bz,double &bphi) const;

    KOKKOS_INLINE_FUNCTION void bmag_interpol(const SimdVector& v, Simd<double>& bmag) const;

    KOKKOS_INLINE_FUNCTION void bmag_interpol(const SimdVector2D& x, const Simd<double>& phi,Simd<double>& bmag) const;

    KOKKOS_INLINE_FUNCTION void follow_field(const SimdVector2D &x_org,const Simd<double>& phi_org, const Simd<double>& phi_dest,SimdVector2D &x_dest) const;
};

#include "magnetic_field.tpp"
#endif