XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
magnetic_field.hpp
Go to the documentation of this file.
1 #ifndef MAGNETIC_FIELD_HPP
2 #define MAGNETIC_FIELD_HPP
3 #include "analytic_psi.hpp"
4 #include "analytic_I.hpp"
6 #include "equil.hpp"
7 #include "bicub.hpp"
8 #include "cub_interp.hpp"
9 #include "tricub.hpp"
10 
11 // Magnetic field class
12 template<class Device>
14  // This line enables private access between the host and device implementations of this class
15  // i.e. Makes the copy from host to device more straightforward
16  template<class Device2> friend class MagneticField;
17 
18  KOKKOS_INLINE_FUNCTION void derivs(const double (&x)[2],double phi,double (&dx)[2]) const;
19 
20 #ifndef STELLARATOR_TRICUB
21  KOKKOS_INLINE_FUNCTION double I_value(double psi_in,int rgn3) const;
22  KOKKOS_INLINE_FUNCTION double I_deriv(double psi_in,int rgn3) const;
23 #endif
24 
25  int ff_step;
26  int ff_order;
27 
28 #ifndef STELLARATOR_TRICUB
31 #else
32  Tricub<Device> psi_tricub;
33  Tricub<Device> Br_tricub;
34  Tricub<Device> Bz_tricub;
35  Tricub<Device> Bphi_tricub;
36 #endif
37 
38  public:
39 
40  double bt_sign;
41  double bp_sign;
43  double inpsi;
44  double outpsi;
46 
47  inline MagneticField(NLReader::NamelistReader& nlr, const MagneticEquilFiles& equil_files, const View<Equil::XPoint*, HostType>& xpts);
48 
49  // Create a mirror with a different device type
50  template<class Device2>
53  m.bt_sign = bt_sign;
54  m.bp_sign = bp_sign;
55  m.ff_step = ff_step;
56  m.ff_order = ff_order;
57  m.bounds = bounds;
58 
59  m.inpsi = inpsi;
60  m.outpsi = outpsi;
61 
62  m.equil = equil;
63 #ifndef STELLARATOR_TRICUB
64  m.psi_bicub = psi_bicub.template mirror<Device2>();
65  m.I_interp = I_interp.template mirror<Device2>();
66 #else
67  m.psi_tricub = psi_tricub.template mirror<Device2>();
68  m.Br_tricub = Br_tricub.template mirror<Device2>();
69  m.Bz_tricub = Bz_tricub.template mirror<Device2>();
70  m.Bphi_tricub = Bphi_tricub.template mirror<Device2>();
71 #endif
72 
73  return m;
74  }
75 
76 #ifndef STELLARATOR_TRICUB
77  /* Constructor for analytic test field
78  * Note safety_factor_coeff is not the safety factor, it just scales with it.
79  * If eq_mz and eq_mpsi are not provided, they will be set to eq_mr.
80  */
81  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)
82  : inpsi(0.0),
83  outpsi(2.0),
84 
85  bt_sign(-1),
86  bp_sign(1),
87 
88  bounds(0.5, 2.5, -1.0, 1.0),
89 
90  ff_step(3),
91  ff_order(4),
92 
93  equil(psi_opt)
94  {
95  // eq_mr = 100 as default
96  int eq_mr = (eq_mr_in==-1 ? 100 : eq_mr_in);
97 
98  /* Set up psi_bicub */
99  // Set eq_mz to eq_mr if not provided
100  int eq_mz = (eq_mz_in==-1 ? eq_mr : eq_mz_in);
101 
102  // Set up linear grid in r and z
103  auto rgrid = create_range_view("rgrid", equil.bounds.min_r, equil.bounds.max_r, eq_mr);
104  auto zgrid = create_range_view("zgrid", equil.bounds.min_z, equil.bounds.max_z, eq_mz);
105 
106  // Get psi and initialize object
107  psi_bicub = Bicub<Device>(rgrid, zgrid, analytic_psi_grid(equil, rgrid, zgrid));
108 
109  /* Set up I_interp */
110  // Set eq_mpsi to eq_mr if not provided
111  int eq_mpsi = (eq_mpsi_in==-1 ? eq_mr : eq_mpsi_in);
112 
113  // Set up linear grid in psi using inpsi and outpsi as psi bounds
114  auto psigrid = create_range_view("psigrid", inpsi, outpsi, eq_mpsi);
115 
116  // Get I and initialize object
117  I_interp = CubInterp<Device>(psigrid, analytic_I_grid(safety_factor_coeff, psigrid));
118  }
119 #endif
120 
121  // Default constructor
123 
124  KOKKOS_INLINE_FUNCTION double psi_norm() const;
125 
126  KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector& v,Simd<double>& psi_out) const;
127 
128  KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector2D& x,const Simd<double>& phi, Simd<double>& psi_out) const;
129 
130  KOKKOS_INLINE_FUNCTION double get_psi(double r, double z, double phi) const;
131 
132  KOKKOS_INLINE_FUNCTION void get_psi_and_derivs(double r, double z, double phi, double& psi, double& dpsidr, double& dpsidz, double& dpsidphi) const;
133 
134  KOKKOS_INLINE_FUNCTION void get_psi_unit_vec(const SimdVector2D& x, double phi, Simd<double>& cosa, Simd<double>& sina) const;
135 
136  KOKKOS_INLINE_FUNCTION double geometry_r(double r) const;
137 
138  KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r,double z,double psi) const;
139  KOKKOS_INLINE_FUNCTION bool is_in_region_1(double r,double z,double psi) const;
140  KOKKOS_INLINE_FUNCTION void check_boundaries(const SimdVector2D& x, Simd<bool>& rz_outside) const;
141  KOKKOS_INLINE_FUNCTION void get_theta(const SimdVector2D& x, Simd<double>& theta) const;
142 
143  // Get complete B-field info
144  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;
145 
146  KOKKOS_INLINE_FUNCTION void bvec_interpol(double r,double z,double phi,double &br,double &bz,double &bphi) const;
147 
148  KOKKOS_INLINE_FUNCTION void bmag_interpol(const SimdVector& v, Simd<double>& bmag) const;
149 
150  KOKKOS_INLINE_FUNCTION void bmag_interpol(const SimdVector2D& x, const Simd<double>& phi,Simd<double>& bmag) const;
151 
152  KOKKOS_INLINE_FUNCTION void follow_field(const SimdVector2D &x_org,const Simd<double>& phi_org, const Simd<double>& phi_dest,SimdVector2D &x_dest) const;
153 };
154 
155 #include "magnetic_field.tpp"
156 #endif
Definition: tricub.hpp:13
double inpsi
Boundary condition used in a few spots.
Definition: magnetic_field.hpp:43
KOKKOS_INLINE_FUNCTION bool is_in_region_1(double r, double z, double psi) const
Definition: magnetic_field.tpp:73
CubInterp< Device > I_interp
The object for interpolating I (for deriving toroidal magnetic field)
Definition: magnetic_field.hpp:30
KOKKOS_INLINE_FUNCTION void get_psi_and_derivs(double r, double z, double phi, double &psi, double &dpsidr, double &dpsidz, double &dpsidphi) const
Definition: magnetic_field.tpp:121
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:321
Definition: simd.hpp:149
Definition: magnetic_equil_files.hpp:16
PsiOption
Definition: equil.hpp:11
KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector &v, Simd< double > &psi_out) const
Definition: magnetic_field.tpp:105
KOKKOS_INLINE_FUNCTION void follow_field(const SimdVector2D &x_org, const Simd< double > &phi_org, const Simd< double > &phi_dest, SimdVector2D &x_dest) const
Definition: magnetic_field.tpp:403
KOKKOS_INLINE_FUNCTION void check_boundaries(const SimdVector2D &x, Simd< bool > &rz_outside) const
Definition: magnetic_field.tpp:78
KOKKOS_INLINE_FUNCTION void get_theta(const SimdVector2D &x, Simd< double > &theta) const
Definition: magnetic_field.tpp:83
KOKKOS_INLINE_FUNCTION double geometry_r(double r) const
Definition: magnetic_field.tpp:183
Definition: rz_bounds.hpp:4
int ff_order
Order of RK scheme used for field following. Can be 1, 2, or 4.
Definition: magnetic_field.hpp:26
View< double *, HostType > create_range_view(std::string name, double x_min, double x_max, int n)
Definition: range_view.hpp:6
Definition: NamelistReader.hpp:193
Definition: magnetic_field.hpp:13
KOKKOS_INLINE_FUNCTION double I_deriv(double psi_in, int rgn3) const
Definition: magnetic_field.tpp:167
RZBounds bounds
Simulation boundary.
Definition: magnetic_field.hpp:42
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)
Definition: magnetic_field.hpp:81
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:45
KOKKOS_INLINE_FUNCTION void bmag_interpol(const SimdVector &v, Simd< double > &bmag) const
Definition: magnetic_field.tpp:372
Definition: cub_interp.hpp:12
MagneticField()
Definition: magnetic_field.hpp:122
KOKKOS_INLINE_FUNCTION double I_value(double psi_in, int rgn3) const
Definition: magnetic_field.tpp:156
KOKKOS_INLINE_FUNCTION double psi_norm() const
Definition: magnetic_field.tpp:62
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
Definition: magnetic_field.tpp:197
double outpsi
Boundary condition used in a few spots.
Definition: magnetic_field.hpp:44
double max_z
Definition: rz_bounds.hpp:8
double bt_sign
Whether toroidal field is reversed?
Definition: magnetic_field.hpp:40
double min_r
Definition: rz_bounds.hpp:5
Definition: bicub.hpp:14
RZBounds bounds
Min and max for r and z.
Definition: equil.hpp:77
double bp_sign
Whether poloidal field is reversed?
Definition: magnetic_field.hpp:41
Definition: simd.hpp:139
KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r, double z, double psi) const
Definition: magnetic_field.tpp:68
double min_z
Definition: rz_bounds.hpp:7
Bicub< Device > psi_bicub
The object for interpolating psi (magnetic flux surfaces)
Definition: magnetic_field.hpp:29
int ff_step
Number of steps taken when projecting the particle location onto the midplane.
Definition: magnetic_field.hpp:25
KOKKOS_INLINE_FUNCTION void get_psi_unit_vec(const SimdVector2D &x, double phi, Simd< double > &cosa, Simd< double > &sina) const
Definition: magnetic_field.tpp:140
KOKKOS_INLINE_FUNCTION void derivs(const double(&x)[2], double phi, double(&dx)[2]) const
Definition: magnetic_field.tpp:517
Definition: equil.hpp:29
MagneticField< Device2 > mirror() const
Definition: magnetic_field.hpp:51
double max_r
Definition: rz_bounds.hpp:6