XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
equil.hpp
Go to the documentation of this file.
1 #ifndef EQUIL_HPP
2 #define EQUIL_HPP
3 #include "space_settings.hpp"
4 #include "simd.hpp"
5 #include "NamelistReader.hpp"
6 #include "rz_bounds.hpp"
7 #include "grid_structs.hpp"
10 
11 // Options for generating analytic equilibrium for testing
12 enum PsiOption{
13  NoXPoints = 0,
16 };
17 
18 namespace Equil{
19 class XPoint{
20  public:
21 
22  double r;
23  double z;
24  double psi;
25  double slope;
26 };
27 }
28 
29 
30 class Equilibrium {
31 
32  public:
33 
34  Equilibrium(NLReader::NamelistReader& nlr, const MagneticEquilFiles::Ptr& equil_files, const View<Equil::XPoint*, HostType>& xpts, const RZPair& axis_in);
35 
36  // Constructor for tests
38  : bounds(0.5, 2.5, -1.0, 1.0),
39 
40  // Axis
41  axis{1.5, 0.1},
42 
43  epsil_psi(1.0e-05),
44 
45  // Xpt1
46  // XGC simulates circular case by putting the xpoint outside the model
47  // Place xpoint below min_z if no X-points
48  // Circular case still uses xpt_psi as a reference, but it's far away; so increase it
49  xpt{1.6, (psi_opt == NoXPoints) ? (bounds.min_z - 1.0) : -0.8},
50  xpt_slope(-0.5),
51  xpt_psi((psi_opt == NoXPoints) ? 5.0 : .8),
53 
54  // Xpt2
55  xpt2{1.4, 0.7},
56  xpt2_slope(-0.5),
57  xpt2_psi(1.0),
58  set_xpt2(psi_opt == TwoXPoints),
59 
60  // Used by profile testing
61  out_decay_factor(0.25),
62  out_decay_width(2.0),
65  {}
66 
67  // Default constructor
69 
70  KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r,double z,double psi) const;
71  KOKKOS_INLINE_FUNCTION bool is_in_region_1(double r,double z,double psi) const;
72  KOKKOS_INLINE_FUNCTION bool is_in_region_3b(double r, double z, double psi) const;
73  KOKKOS_INLINE_FUNCTION void check_boundaries(const SimdVector2D& x, Simd<bool>& rz_outside) const;
74  KOKKOS_INLINE_FUNCTION void get_theta(const SimdVector2D& x, Simd<double>& theta) const;
75 
76  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,
77  const View<double**,CLayout,HostType>& eq_B_R_rz = View<double**,CLayout,HostType>(NoInit("eq_B_R_rz"), 0, 0),
78  const View<double**,CLayout,HostType>& eq_B_Z_rz = View<double**,CLayout,HostType>(NoInit("eq_B_Z_rz"), 0, 0),
79  const View<double**,CLayout,HostType>& eq_B_phi_rz = View<double**,CLayout,HostType>(NoInit("eq_B_phi_rz"), 0, 0),
80  const View<double**,CLayout,HostType>& eq_VMEC_angle_rz = View<double**,CLayout,HostType>(NoInit("eq_VMEC_angle_rz"), 0, 0),
81  const View<double**,CLayout,HostType>& eq_PEST_angle_rz = View<double**,CLayout,HostType>(NoInit("eq_PEST_angle_rz"), 0, 0));
82 
83  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);
84 
85  // Equil variables
87  double xpt_psi;
88  double psi_norm;
89  double epsil_psi = 1.0e-05;
91  double xpt_slope;
92  bool set_xpt2;
94  double xpt2_slope;
95  double xpt2_psi;
97  double axis_b;
98 
99  // Use in eq profile calculations
104 };
105 
106 #include "equil.tpp"
107 
108 #endif
Definition: equil.hpp:14
PsiOption
Definition: equil.hpp:12
KOKKOS_INLINE_FUNCTION bool is_in_region_1(double r, double z, double psi) const
Definition: equil.tpp:29
Definition: equil.hpp:19
KOKKOS_INLINE_FUNCTION bool is_in_region_3b(double r, double z, double psi) const
Definition: equil.tpp:45
Equilibrium()
Definition: equil.hpp:68
Definition: rz_bounds.hpp:4
RZPair xpt2
coordinates of 2nd X-point
Definition: equil.hpp:93
Definition: NamelistReader.hpp:193
double z
Definition: equil.hpp:23
double out_decay_width
width for exponential decay for psi&gt;sml_outpsi
Definition: equil.hpp:102
double axis_b
magnetic field magnitude on axis
Definition: equil.hpp:97
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)
Definition: equil.cpp:130
Definition: grid_structs.hpp:28
bool set_xpt2
Whether to use a 2nd X-point.
Definition: equil.hpp:92
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))
Definition: equil.cpp:70
double out_decay_factor
profiles decay exponentially to f(sml_outpsi)/decay_factor for psi&gt;sml_outpsi
Definition: equil.hpp:100
double priv_flux_decay_factor
profiles decay exponentially to f(sml_outpsi)/decay_factor in priv. flux region
Definition: equil.hpp:101
RZBounds bounds
Min and max for r and z.
Definition: equil.hpp:86
double psi_norm
Psi value to use for normalization.
Definition: equil.hpp:88
double xpt2_psi
psi coordinate at 2nd X-point
Definition: equil.hpp:95
double xpt_slope
Slope (which slope?) at 1st X-point.
Definition: equil.hpp:91
double xpt2_slope
Slope (which slope?) at 2nd X-point.
Definition: equil.hpp:94
KOKKOS_INLINE_FUNCTION void check_boundaries(const SimdVector2D &x, Simd< bool > &rz_outside) const
Definition: equil.tpp:56
double psi
Definition: equil.hpp:24
RZPair xpt
coordinates of 1st X-point
Definition: equil.hpp:90
KOKKOS_INLINE_FUNCTION void get_theta(const SimdVector2D &x, Simd< double > &theta) const
Definition: equil.tpp:66
Definition: simd.hpp:18
Definition: equil.hpp:15
Definition: simd.hpp:139
Equilibrium(PsiOption psi_opt)
Definition: equil.hpp:37
double min_z
Definition: rz_bounds.hpp:7
double epsil_psi
Not sure?
Definition: equil.hpp:89
RZPair axis
coordinates of axis
Definition: equil.hpp:96
double r
Definition: equil.hpp:22
Definition: equil.hpp:13
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:69
double slope
Definition: equil.hpp:25
Definition: equil.hpp:30
KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r, double z, double psi) const
Definition: equil.tpp:6
std::shared_ptr< MagneticEquilFiles > Ptr
Definition: magnetic_equil_files.hpp:11
double priv_flux_decay_width
width for exponential decay in private flux region
Definition: equil.hpp:103
double xpt_psi
Psi coordinate of 1st X-point.
Definition: equil.hpp:87