magnetic_field.hpp

Go to the documentation of this file.

#ifndef MAGNETIC_FIELD_HPP
#define MAGNETIC_FIELD_HPP
#include "analytic_psi.hpp"
#include "analytic_I.hpp"
#include "magnetic_equil_files.hpp"
#include "equil.hpp"
#include "bicub.hpp"
#include "cub_interp.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_TRICUB
    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_TRICUB
    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; 

    inline MagneticField(NLReader::NamelistReader& nlr, const MagneticEquilFiles& equil_files, const View<Equil::XPoint*, HostType>& xpts);

    // 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_TRICUB
        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;
    }

#ifndef STELLARATOR_TRICUB
    /* Constructor for analytic test field
     * Note safety_factor_coeff is not the safety factor, it just scales with it.
     * If eq_mz and eq_mpsi are not provided, they will be set to eq_mr.
     */
    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)
      : inpsi(0.0),
        outpsi(2.0),

        bt_sign(-1),
        bp_sign(1),

        bounds(0.5, 2.5, -1.0, 1.0),

        ff_step(3),
        ff_order(4),

        equil(psi_opt)
    {
        // eq_mr = 100 as default
        int eq_mr = (eq_mr_in==-1 ? 100 : eq_mr_in);

        /* Set up psi_bicub */
        // Set eq_mz to eq_mr if not provided
        int eq_mz = (eq_mz_in==-1 ? eq_mr : eq_mz_in);

        // Set up linear grid in r and z
        auto rgrid = create_range_view("rgrid", equil.bounds.min_r, equil.bounds.max_r, eq_mr);
        auto zgrid = create_range_view("zgrid", equil.bounds.min_z, equil.bounds.max_z, eq_mz);

        // Get psi and initialize object
        psi_bicub = Bicub<Device>(rgrid, zgrid, analytic_psi_grid(equil, rgrid, zgrid));

        /* Set up I_interp */
        // Set eq_mpsi to eq_mr if not provided
        int eq_mpsi = (eq_mpsi_in==-1 ? eq_mr : eq_mpsi_in);

        // Set up linear grid in psi using inpsi and outpsi as psi bounds
        auto psigrid = create_range_view("psigrid", inpsi, outpsi, eq_mpsi);

        // Get I and initialize object
        I_interp = CubInterp<Device>(psigrid, analytic_I_grid(safety_factor_coeff, psigrid));
    }
#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_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;

    // Get complete B-field info
    KOKKOS_INLINE_FUNCTION void field(const SimdVector& v, SimdVector &bvec, SimdVector (&jacb)[3], Simd<double>& psivec, SimdVector2D &gradpsi, SimdVector &tdb, 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