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"

// Magnetic field class
template<class Device>
class MagneticField {

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

    public:

    // Field variables
    double bt_sign;  
    double bp_sign;  
    int ff_step;     
    int ff_order;    
    RZBounds bounds; 

    double inpsi;    
    double outpsi;   


    // Structures inside B-field
    Bicub<Device> psi_bicub;   
    CubInterp<Device> I_interp;   
    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;

        m.psi_bicub = psi_bicub.template mirror<Device2>();
        m.I_interp = I_interp.template mirror<Device2>();

        return m;
    }

    /* 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));
    }

    // Default constructor
    MagneticField(){}

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

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

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

    KOKKOS_INLINE_FUNCTION double I_value(double psi_in,int rgn3) const;

    KOKKOS_INLINE_FUNCTION double I_deriv(double psi_in,int rgn3) const;

    KOKKOS_INLINE_FUNCTION double geometry_r(double r) const;

    // Get complete B-field info
    KOKKOS_INLINE_FUNCTION void field(const SimdVector2D& x, 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