XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
grid_field_pack.hpp
Go to the documentation of this file.
1 #ifndef GRID_FIELD_PACK_HPP
2 #define GRID_FIELD_PACK_HPP
3 #include "space_settings.hpp"
4 #include "NamelistReader.hpp"
5 
6 #include "globals.hpp"
7 #include "grid.hpp"
8 #include "linear_weights.hpp"
9 #include "local_fields.hpp"
10 #include "push_controls.hpp"
11 #include "field_weights.hpp"
12 #include "gyro_radius.hpp"
13 #include "grid_field.hpp"
14 
15 // Base struct for GridFieldPack so that non-GPU code doesn't need to be templated just because there is a GridFieldPack argument
17 };
18 
19 template<class Device, PhiInterpType PIT>
21  private:
24 
25  public:
26 
27  bool turb_efield; // Ideally this would not be here but it is tedious to pass into the gather otherwise
28 
29  int phi_offset;
31 
32  int nrho;
33  double inv_drho;
34 
35  // Device Views
36 #ifdef XGC1
40 # ifdef DELTAF_CONV
43 # endif
44 # ifdef EXPLICIT_EM
45  //For EM
56 # endif
57 #else
60 # ifdef SONIC_GK
64 # endif
65 #endif
66 
68  : is_initialized(false) {}
69 
70  GridFieldPack(KinType kintype_in, bool turb_efield_in, int nrho_in, double inv_drho_in)
71  : is_initialized(true),
72  kintype(kintype_in),
73  turb_efield(turb_efield_in),
74  nrho(nrho_in),
75  inv_drho(inv_drho_in),
76  phi_offset(0),
77  node_offset(0) {}
78 
79  bool can_reuse(KinType kintype_in){
80  if(is_initialized){
81 #ifdef XGC1
82  // For XGC1, can reuse if the existing GridFieldPack is the same KinType
83  return kintype_in==kintype;
84 #else
85  // For XGCa, the drift kinetic GridFieldPack is a subset of the gyrokinetic one,
86  // so is always reused
87  return true;
88 #endif
89  }else{
90  // Cant reuse if uninitialized
91  return false;
92  }
93  }
94 
95  // Device functions:
96 
97  KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector& vec, int i_simd, const SimdVector& B, double Bmag) const{
98  vec.phi[i_simd]=(vec.phi[i_simd]*Bmag - vec.r[i_simd]*B.r[i_simd] - vec.z[i_simd]*B.z[i_simd]) / B.phi[i_simd];
99  }
100 
101 
102 
103  KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int i_node, double wp, const FieldCorrection& corr, const FieldWeights<DriftKin, PIT>& wts, LocalFields& fld) const {
104  const int node = i_node - node_offset;
105 #ifdef XGCA
106  E_rho(node,0).gather(fld.E,i_simd, wp, corr); // irho = 0
107 #else
108  const int i_phi = (wts.phi.i - phi_offset)%E_phi_ff.f.extent(0); //should be grid.nplanes (?);
109  E_phi_ff(i_phi,node).gather(fld.E,i_simd, wp, wts.phi.w, corr);
110 #ifdef EXPLICIT_EM
111  Ah_phi_ff(i_phi,node).gather(fld.Ah,i_simd, wp, wts.phi.w);
112  dAh_phi_ff(i_phi,node).gather(fld.dAh,i_simd, wp, wts.phi.w, corr);
113 
114  if(push_controls.em_mixed_variable) {
115  dAs_phi_ff(i_phi,node).gather(fld.dAs,i_simd, wp, wts.phi.w, corr);
116  As_phi_ff(i_phi,node).gather(fld.As,i_simd, wp, wts.phi.w);
117  }
118  if(push_controls.em_mixed_variable && (push_controls.em_pullback_mode==4)) {
119  Epar_em_phi_ff(i_phi,node).gather(fld.Epar_em,i_simd,wp,wts.phi.w);
120  }
121 #endif
122 #endif
123 
124 #if defined(XGC1) && defined(DELTAF_CONV)
125  E00_ff(node).gather(fld.E00,i_simd, wp, wts.phi.w, corr);
126  ddpotdt_phi_ff(i_phi,node).gather(fld.ddpotdt,i_simd, wp, wts.phi.w);
127 #endif
128  }
129 
130  // Gather all fields (ion)
131  KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int node, double wp, const FieldCorrection& corr, const FieldWeights<GyroKin, PIT>& wts, LocalFields& fld) const {
132 #ifdef NEWGYROMATRIX
133  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
134  int irho = wts.rho[0].i;
135 #else
136  int irho = wts.rho.i;
137 #endif
138 #ifdef XGCA
139  E_rho(node,irho).gather(fld.E,i_simd, wp, wts, corr, E_rho(node,irho+1));
140 # ifdef SONIC_GK
141  dEr_B2_rho(node,0).gather(fld.dEr_B2,i_simd, wp, corr); // irho = 0
142  dEz_B2_rho(node,0).gather(fld.dEz_B2,i_simd, wp, corr); // irho = 0
143  du2_E_rho(node,0).gather(fld.du2_E,i_simd, wp, corr); // irho = 0
144 # endif
145 #else
146  E_rho_ff(node,irho).gather(fld.E,i_simd, wp, wts, corr, E_rho_ff(node,irho+1));
147 #ifdef EXPLICIT_EM
148  Ah_rho_ff(node,irho).gather(fld.Ah,i_simd, wp, wts, Ah_rho_ff(node,irho+1));
149  dAh_rho_ff(node,irho).gather(fld.dAh,i_simd, wp, wts, corr, dAh_rho_ff(node,irho+1));
150 
151  if(push_controls.em_mixed_variable) {
152  dAs_rho_ff(node,irho).gather(fld.dAs,i_simd, wp, wts, corr, dAs_rho_ff(node,irho+1));
153  As_rho_ff(node,irho).gather(fld.As,i_simd, wp, wts, As_rho_ff(node,irho+1));
154  }
155  if(push_controls.em_mixed_variable && (push_controls.em_pullback_mode==4)) {
156  Epar_em_rho_ff(node,irho).gather(fld.Epar_em,i_simd,wp,wts, Epar_em_rho_ff(node,irho+1));
157  }
158 #endif
159 #endif
160  }
161 
162 
174  template<KinType PT>
175  KOKKOS_INLINE_FUNCTION void fields_at_point(const PushControls &push_controls, const Grid<Device> &grid, const SimdVector &B, const SimdVector2D &gradpsi, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<PT>& rho, LocalFields& fld) const
176  {
177  // Initialize local fields to 0
178  fld.E.zero();
179 #ifdef DELTAF_CONV
180  fld.E00.zero();
181  fld.ddpotdt.zero();
182 #endif
183 #ifdef EXPLICIT_EM
184  fld.dAh.zero();
185  fld.dAs.zero();
186  fld.Ah.zero();
187  fld.As.zero();
188  fld.Epar_em.zero();
189 #endif
190 #ifdef SONIC_GK
191  fld.dEr_B2.zero();
192  fld.dEz_B2.zero();
193  fld.du2_E.zero();
194 #endif
195 
196  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
197  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
198  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
199 
200  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
201 
202  // Get rho and phi weights
203  FieldWeights<PT, PIT> wts(grid, phi, i_simd, rho, inv_drho, nrho);
204 
205  // Field basis correction
206  FieldCorrection field_correction(gradpsi, i_simd);
207 
208  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
209  for (int ip = 0; ip<3; ip++){
210  int node = itr[i_simd]>0 ? grid.get_node_index(itr_work, ip) : 0;
211  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
212 
213  field_correction.set(i_simd, grid.uses_rz_basis(node), gradpsi);
214 
215  gather_all_fields(push_controls, i_simd, node, wp, field_correction, wts, fld);
216  }
217 
218 #ifdef NEOCLASSICAL_TEST
219  fld.E.phi[i_simd]=0.0;
220 #else
221 #ifdef XGC1
222  if(turb_efield){
223  // Correct phi component of E-field
224  double Bmag = B.magnitude(i_simd);
225 
226  // Get phi from para
227  phi_from_para(fld.E, i_simd, B, Bmag);
228 #ifdef EXPLICIT_EM
229  phi_from_para(fld.dAh, i_simd, B, Bmag);
230  phi_from_para(fld.dAs, i_simd, B, Bmag);
231 #endif
232  }else{
233  fld.E.phi[i_simd]=0.0;
234 #ifdef EXPLICIT_EM
235  fld.dAh.phi[i_simd]=0.0;
236  fld.dAs.phi[i_simd]=0.0;
237 #endif
238  }
239 #endif
240 #endif
241  }
242  }
243 
244 
245  KOKKOS_INLINE_FUNCTION void get_Ah(const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<GyroKin>& rho, Simd<double>& Ah) const
246  {
247  // Initialize local fields to 0
248 #ifdef EXPLICIT_EM
249  Ah = 0.0;
250 
251  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
252  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
253  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
254  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
255 
256  // Get rho and phi weights
257  FieldWeights<GyroKin, PIT> wts(grid, phi, i_simd, rho, inv_drho, nrho);
258 #ifdef NEWGYROMATRIX
259  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
260  int irho = wts.rho[0].i;
261 #else
262  int irho = wts.rho.i;
263 #endif
264 
265  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
266  for (int ip = 0; ip<3; ip++){
267  int node = itr[i_simd]>0 ? grid.get_node_index(itr_work, ip) : 0;
268  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
269 
270  Ah_rho_ff(node,irho).gather(Ah,i_simd, wp, wts, Ah_rho_ff(node,irho+1));
271  }
272  }
273 #endif
274  }
275 
276  KOKKOS_INLINE_FUNCTION void get_Ah(const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<DriftKin>& rho, Simd<double>& Ah) const
277  {
278  // Initialize local fields to 0
279 #ifdef EXPLICIT_EM
280  Ah = 0.0;
281 
282  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
283  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
284  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
285  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
286 
287  // Get rho and phi weights
288  FieldWeights<DriftKin, PIT> wts(grid, phi, i_simd, rho, inv_drho, nrho);
289 
290  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
291  for (int ip = 0; ip<3; ip++){
292  int i_node = itr[i_simd]>0 ? grid.get_node_index(itr_work, ip) : 0;
293  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
294 
295  const int i_phi = (wts.phi.i - phi_offset)%grid.nplanes;
296  const int node = i_node - node_offset;
297 
298  Ah_phi_ff(i_phi,node).gather(Ah,i_simd, wp, wts.phi.w);
299  }
300  }
301 #endif
302  }
303 
304  KOKKOS_INLINE_FUNCTION void get_dpot(const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, Simd<double>& dpot_out) const {
305 
306  dpot_out.zero();
307 
308  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
309  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
310  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
311 
312  double phi = itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
314  FieldWeights<DriftKin, PIT> wts(grid, phi, i_simd, rho, inv_drho, nrho);
315 
316  for(int ip=0; ip<3; ip++){
317  int node = itr[i_simd]>0 ? grid.get_node_index(itr_work, ip) : 0;
318  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
319 
320 #ifdef XGC1
321  //#ifndef EXPLICIT_EM
322  // for electron -- rho=0 case, optimized.
323  dpot_ff(node).gather(dpot_out,i_simd, wp, wts.phi.w);
324  //#else
325  // dpot_out[i_simd]+=wp*dpot_n0(node);
326  //#endif
327 #elif defined(XGCA)
328  // Simple gather
329  dpot(node).gather(dpot_out,i_simd,wp);
330 #endif
331  }
332  }
333  }
334 };
335 #endif
KOKKOS_INLINE_FUNCTION void set(int i_simd, double basis, const SimdVector2D &gradpsi)
Definition: field.hpp:40
Simd< double > r
Definition: simd.hpp:150
Definition: simd.hpp:149
GridFieldPack()
Definition: grid_field_pack.hpp:67
Definition: gyro_radius.hpp:24
KOKKOS_INLINE_FUNCTION void get_dpot(const Grid< Device > &grid, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, Simd< double > &dpot_out) const
Definition: grid_field_pack.hpp:304
KOKKOS_INLINE_FUNCTION void get_Ah(const Grid< Device > &grid, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, SimdGyroRadius< GyroKin > &rho, Simd< double > &Ah) const
Definition: grid_field_pack.hpp:245
bool turb_efield
Definition: grid_field_pack.hpp:27
double inv_drho
1.0 over rho spacing - Used if kintype==Gyrokin
Definition: grid_field_pack.hpp:33
Definition: field_weights.hpp:13
Definition: globals.hpp:82
Definition: push_controls.hpp:8
int phi_offset
Offset for phi_ff field decomposition.
Definition: grid_field_pack.hpp:29
GridField< Device, vec2d_if_axisym< PIT >), PIT, TorType::OnePlane, KinType::GyroKin > E_rho
Definition: grid_field_pack.hpp:58
Definition: grid.hpp:53
int nplanes
Number of planes.
Definition: grid.hpp:217
Definition: grid_field_pack.hpp:20
KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int i_node, double wp, const FieldCorrection &corr, const FieldWeights< DriftKin, PIT > &wts, LocalFields &fld) const
Definition: grid_field_pack.hpp:103
KOKKOS_INLINE_FUNCTION int uses_rz_basis(const int inode) const
Definition: grid.tpp:900
Definition: local_fields.hpp:7
Definition: grid_field.hpp:25
bool can_reuse(KinType kintype_in)
Definition: grid_field_pack.hpp:79
GridFieldPack(KinType kintype_in, bool turb_efield_in, int nrho_in, double inv_drho_in)
Definition: grid_field_pack.hpp:70
KOKKOS_INLINE_FUNCTION void get_Ah(const Grid< Device > &grid, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, SimdGyroRadius< DriftKin > &rho, Simd< double > &Ah) const
Definition: grid_field_pack.hpp:276
KinType kintype
Definition: grid_field_pack.hpp:22
int nrho
Gyroavg dimension (not including 0) - Used if kintype==GyroKin.
Definition: grid_field_pack.hpp:32
KOKKOS_INLINE_FUNCTION int get_node_index(int triangle_index, int tri_vertex_index) const
Definition: grid.tpp:777
Definition: field.hpp:28
KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector &vec, int i_simd, const SimdVector &B, double Bmag) const
Definition: grid_field_pack.hpp:97
Definition: globals.hpp:83
Simd< double > phi
Definition: simd.hpp:152
int node_offset
Offset for phi_ff field decomposition.
Definition: grid_field_pack.hpp:30
Simd< double > z
Definition: simd.hpp:151
KinType
Definition: globals.hpp:81
KOKKOS_INLINE_FUNCTION void zero()
Definition: simd.hpp:76
Definition: grid_structs.hpp:7
GridField< Device, VarType::Scalar, PIT, TorType::OnePlane, KinType::DriftKin > dpot
Definition: grid_field_pack.hpp:59
Definition: simd.hpp:139
KOKKOS_INLINE_FUNCTION void zero()
Definition: simd.hpp:167
bool is_initialized
Definition: grid_field_pack.hpp:23
Definition: gyro_radius.hpp:96
KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int node, double wp, const FieldCorrection &corr, const FieldWeights< GyroKin, PIT > &wts, LocalFields &fld) const
Definition: grid_field_pack.hpp:131
Definition: gyro_radius.hpp:74
KOKKOS_INLINE_FUNCTION double magnitude(const int i_simd) const
Definition: simd.hpp:176
Definition: grid_field_pack.hpp:16
SimdVector E
Definition: local_fields.hpp:8
KOKKOS_INLINE_FUNCTION void fields_at_point(const PushControls &push_controls, const Grid< Device > &grid, const SimdVector &B, const SimdVector2D &gradpsi, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, SimdGyroRadius< PT > &rho, LocalFields &fld) const
Definition: grid_field_pack.hpp:175