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 "magnetic_equil_files.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:

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

    inline Equilibrium(double eq_min_r, double eq_max_r, double eq_min_z, double eq_max_z,
                       double eq_x_psi, double epsil_psi_in, double eq_x_r, double eq_x_slope, double eq_x_z,
                       double eq_x2_r, double eq_x2_slope, double eq_x2_z, double eq_x2_psi, double eq_axis_r,
                       double eq_axis_z, int eq_mpsi);

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

        // Axis
        axis_r(1.5),
        axis_z(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_r(1.6),
        xpt_z((psi_opt == NoXPoints) ? (bounds.min_z - 1.0) : -0.8),
        xpt_slope(-0.5),
        xpt_psi((psi_opt == NoXPoints) ? 5.0 : .8),

        // Xpt2
        xpt2_r(1.4),
        xpt2_z(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 void check_boundaries(const SimdVector2D& x, Simd<bool>& rz_outside) const;
    KOKKOS_INLINE_FUNCTION void get_theta(const SimdVector2D& x, Simd<double>& theta) const;

    inline 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 epsil_psi = 1.0e-05;  
    double xpt_r;      
    double xpt_slope;  
    double xpt_z;      
    bool set_xpt2;     
    double xpt2_r;     
    double xpt2_slope; 
    double xpt2_z;     
    double xpt2_psi;   
    double axis_r;     
    double axis_z;     

    // 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