equil.hpp

Go to the documentation of this file.

#ifndef EQUIL_HPP
#define EQUIL_HPP
#include "space_settings.hpp"
#include "simd.hpp"
#include "NamelistReader.hpp"
#include "rz_bounds.hpp"
#include "grid_structs.hpp"
#include "magnetic_equil_files.hpp"
#include "domain_decomposition.hpp"

// Options for generating analytic equilibrium for testing
enum PsiOption{
    NoXPoints = 0,
    OneXPoint,
    TwoXPoints
};

namespace Equil{
class XPoint{
    public:

    double r;
    double z;
    double psi;
    double slope;
};
}


class Equilibrium {

    public:

    Equilibrium(NLReader::NamelistReader& nlr, const MagneticEquilFiles::Ptr& equil_files, const View<Equil::XPoint*, HostType>& xpts, const RZPair& axis_in);

    // Constructor for tests
    Equilibrium(PsiOption psi_opt)
      : bounds(0.5, 2.5, -1.0, 1.0),

        // Axis
        axis{1.5, 0.1},

        epsil_psi(1.0e-05),

        // Xpt1
        // XGC simulates circular case by putting the xpoint outside the model
        // Place xpoint below min_z if no X-points
        // Circular case still uses xpt_psi as a reference, but it's far away; so increase it
        xpt{1.6, (psi_opt == NoXPoints) ? (bounds.min_z - 1.0) : -0.8},
        xpt_slope(-0.5),
        xpt_psi((psi_opt == NoXPoints) ? 5.0 : .8),
        psi_norm(xpt_psi),

        // Xpt2
        xpt2{1.4, 0.7},
        xpt2_slope(-0.5),
        xpt2_psi(1.0),
        set_xpt2(psi_opt == TwoXPoints),

        // Used by profile testing
        out_decay_factor(0.25),
        out_decay_width(2.0),
        priv_flux_decay_factor(0.3),
        priv_flux_decay_width(2.5)
    {}

    // Default constructor
    Equilibrium(){}

    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 bool is_in_region_3b(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;

    void write(const DomainDecomposition<DeviceType>& pol_decomp, const View<double*,CLayout,HostType>& eq_psi_grid, const View<double*,CLayout,HostType>& eq_I, const View<double**,CLayout,HostType>& eq_psi_rz, double bt_sign, double bp_sign, bool is_3D, int periodicity = 1, int eq_mp = 1,
               const View<double**,CLayout,HostType>& eq_B_R_rz = View<double**,CLayout,HostType>(NoInit("eq_B_R_rz"), 0, 0),
               const View<double**,CLayout,HostType>& eq_B_Z_rz = View<double**,CLayout,HostType>(NoInit("eq_B_Z_rz"), 0, 0),
               const View<double**,CLayout,HostType>& eq_B_phi_rz = View<double**,CLayout,HostType>(NoInit("eq_B_phi_rz"), 0, 0),
               const View<double**,CLayout,HostType>& eq_VMEC_angle_rz = View<double**,CLayout,HostType>(NoInit("eq_VMEC_angle_rz"), 0, 0),
               const View<double**,CLayout,HostType>& eq_PEST_angle_rz = View<double**,CLayout,HostType>(NoInit("eq_PEST_angle_rz"), 0, 0));

    void set_decay_factors(double out_decay_factor_in, double priv_flux_decay_factor_in, double out_decay_width_in, double priv_flux_decay_width_in, bool set_xpt2_in);

    // Equil variables
    RZBounds bounds;   
    double xpt_psi;    
    double psi_norm;   
    double epsil_psi = 1.0e-05;  
    RZPair xpt;        
    double xpt_slope;  
    bool set_xpt2;     
    RZPair xpt2;       
    double xpt2_slope; 
    double xpt2_psi;   
    RZPair axis;       
    double axis_b;     

    // Use in eq profile calculations
    double out_decay_factor;       
    double priv_flux_decay_factor; 
    double out_decay_width;        
    double priv_flux_decay_width;  
};

#include "equil.tpp"

#endif