XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
14 /* VarType specifies whether the field is a vector, a 2D vector, or a scalar
15  * */
16 enum class VarType{
17  Vector,
18  Vector2D,
19  Scalar
20 };
21 
22 template<VarType VT, PhiInterpType PIT>
23 struct Field{};
24 
25 //Struct used for E_phi_ff, dAs_phi_ff, dAh_phi_ff, As_phi_ff, Ah_phi_ff. 3 vector components, 2 contributing planes
26 template<>
28  double V[3][2];
29 
38  KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const double (&wphi)[2], const double (&rvec)[2], const double (&zvec)[2]) const{
39  vec.r[i_simd] += wp*(wphi[0]*( rvec[0]*V[PIR][0] + zvec[0]*V[PIZ][0] )
40  + wphi[1]*( rvec[0]*V[PIR][1] + zvec[0]*V[PIZ][1] ) );
41  vec.z[i_simd] += wp*(wphi[0]*( rvec[1]*V[PIR][0] + zvec[1]*V[PIZ][0] )
42  + wphi[1]*( rvec[1]*V[PIR][1] + zvec[1]*V[PIZ][1] ) );
43  vec.phi[i_simd] += wp*(wphi[0]* V[PIP][0]
44  + wphi[1]* V[PIP][1] );
45  }
46 
56  KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const FieldWeights<GyroKin, PhiInterpType::Planes>& wts, const double (&rvec)[2], const double (&zvec)[2], const Field<VarType::Vector,PhiInterpType::Planes>& field2) const{
57 # ifdef NEWGYROMATRIX
58  vec.r[i_simd] += wp*(wts.phi.w[0]*( rvec[0]*(wts.rho[0].w[0]*V[PIR][0] + wts.rho[0].w[1]*field2.V[PIR][0])
59  + zvec[0]*(wts.rho[0].w[0]*V[PIZ][0] + wts.rho[0].w[1]*field2.V[PIZ][0]) )
60  + wts.phi.w[1]*( rvec[0]*(wts.rho[1].w[0]*V[PIR][1] + wts.rho[1].w[1]*field2.V[PIR][1])
61  + zvec[0]*(wts.rho[1].w[0]*V[PIZ][1] + wts.rho[1].w[1]*field2.V[PIZ][1]) ) );
62  vec.z[i_simd] += wp*(wts.phi.w[0]*( rvec[1]*(wts.rho[0].w[0]*V[PIR][0] + wts.rho[0].w[1]*field2.V[PIR][0])
63  + zvec[1]*(wts.rho[0].w[0]*V[PIZ][0] + wts.rho[0].w[1]*field2.V[PIZ][0]) )
64  + wts.phi.w[1]*( rvec[1]*(wts.rho[1].w[0]*V[PIR][1] + wts.rho[1].w[1]*field2.V[PIR][1])
65  + zvec[1]*(wts.rho[1].w[0]*V[PIZ][1] + wts.rho[1].w[1]*field2.V[PIZ][1]) ) );
66  vec.phi[i_simd] += wp*(wts.phi.w[0]* (wts.rho[0].w[0]*V[PIP][0] + wts.rho[0].w[1]*field2.V[PIP][0])
67  + wts.phi.w[1]* (wts.rho[1].w[0]*V[PIP][1] + wts.rho[1].w[1]*field2.V[PIP][1]) );
68 # else
69  vec.r[i_simd] += wp*(wts.phi.w[0]*( rvec[0]*(wts.rho.w[0]*V[PIR][0] + wts.rho.w[1]*field2.V[PIR][0])
70  + zvec[0]*(wts.rho.w[0]*V[PIZ][0] + wts.rho.w[1]*field2.V[PIZ][0]) )
71  + wts.phi.w[1]*( rvec[0]*(wts.rho.w[0]*V[PIR][1] + wts.rho.w[1]*field2.V[PIR][1])
72  + zvec[0]*(wts.rho.w[0]*V[PIZ][1] + wts.rho.w[1]*field2.V[PIZ][1]) ) );
73  vec.z[i_simd] += wp*(wts.phi.w[0]*( rvec[1]*(wts.rho.w[0]*V[PIR][0] + wts.rho.w[1]*field2.V[PIR][0])
74  + zvec[1]*(wts.rho.w[0]*V[PIZ][0] + wts.rho.w[1]*field2.V[PIZ][0]) )
75  + wts.phi.w[1]*( rvec[1]*(wts.rho.w[0]*V[PIR][1] + wts.rho.w[1]*field2.V[PIR][1])
76  + zvec[1]*(wts.rho.w[0]*V[PIZ][1] + wts.rho.w[1]*field2.V[PIZ][1]) ) );
77  vec.phi[i_simd] += wp*(wts.phi.w[0]* (wts.rho.w[0]*V[PIP][0] + wts.rho.w[1]*field2.V[PIP][0])
78  + wts.phi.w[1]* (wts.rho.w[0]*V[PIP][1] + wts.rho.w[1]*field2.V[PIP][1]) );
79 # endif
80  }
81 
82 };
83 
84 //Struct used for E00_ff. 2 vector components, 2 contributing planes
85 template<>
87  double V[2][2];
88 
97  KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const double (&wphi)[2], const double (&rvec)[2], const double (&zvec)[2]) const{
98  vec.r[i_simd] += wp*(wphi[0]*( rvec[0]*V[PIR][0] + zvec[0]*V[PIZ][0] )
99  + wphi[1]*( rvec[0]*V[PIR][1] + zvec[0]*V[PIZ][1] ) );
100  vec.z[i_simd] += wp*(wphi[0]*( rvec[1]*V[PIR][0] + zvec[1]*V[PIZ][0] )
101  + wphi[1]*( rvec[1]*V[PIR][1] + zvec[1]*V[PIZ][1] ) );
102  }
103 };
104 
105 //Struct used for ddpotdt_phi_ff, As_phi_ff, Ah_phi_ff, Ah_cv for EM. 2 contributing planes
106 template<>
108  double S[2];
109 
116  KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp, const double (&wphi)[2]) const{
117  sca[i_simd] += wp* ( wphi[0]*S[0] + wphi[1]*S[1] );
118  }
119 
127  KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp, const FieldWeights<GyroKin, PhiInterpType::Planes>& wts, const Field<VarType::Scalar,PhiInterpType::Planes>& field2) const{
128 # ifdef NEWGYROMATRIX
129  sca[i_simd] += wp* ( wts.phi.w[0]*wts.rho[0].w[0]*S[0]
130  + wts.phi.w[1]*wts.rho[1].w[0]*S[1]
131  + wts.phi.w[0]*wts.rho[0].w[1]*field2.S[0]
132  + wts.phi.w[1]*wts.rho[1].w[1]*field2.S[1]);
133 # else
134  sca[i_simd] += wp* ( wts.phi.w[0]*wts.rho.w[0]*S[0]
135  + wts.phi.w[1]*wts.rho.w[0]*S[1]
136  + wts.phi.w[0]*wts.rho.w[1]*field2.S[0]
137  + wts.phi.w[1]*wts.rho.w[1]*field2.S[1]);
138 # endif
139  }
140 };
141 
142 // Struct used for E_rho and du2_E_rho. 2 vector components (r,z)
143 template<>
145  double E[2];
146 
154  KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const double (&rvec)[2], const double (&zvec)[2]) const{
155  vec.r[i_simd] += wp*( rvec[0]*E[PIR] + zvec[0]*E[PIZ] );
156  vec.z[i_simd] += wp*( rvec[1]*E[PIR] + zvec[1]*E[PIZ] );
157  }
158 
168  KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const FieldWeights<GyroKin, PhiInterpType::None>& wts, const double (&rvec)[2], const double (&zvec)[2], const Field<VarType::Vector2D,PhiInterpType::None>& field2) const{
169  vec.r[i_simd] += wp*( rvec[0]*(wts.rho.w[0]*E[PIR] + wts.rho.w[1]*field2.E[PIR])
170  +zvec[0]*(wts.rho.w[0]*E[PIZ] + wts.rho.w[1]*field2.E[PIZ]) );
171  vec.z[i_simd] += wp*( rvec[1]*(wts.rho.w[0]*E[PIR] + wts.rho.w[1]*field2.E[PIR])
172  +zvec[1]*(wts.rho.w[0]*E[PIZ] + wts.rho.w[1]*field2.E[PIZ]) );
173  }
174 
175 };
176 
177 // Struct for scalar without planar interpolation (i.e. just a double)
178 template<>
180  double S;
181 
187  KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp){
188  sca[i_simd] += wp*S;
189  }
190 };
191 
192 // Fortran calls
193 #ifdef XGC1
194 extern "C" void set_rho_ff_pointers(int nrho, int n, double* pot0, double* dpot, Field<VarType::Scalar,PhiInterpType::Planes>* dpot_ff, Field<VarType::Scalar,PhiInterpType::Planes>* pot_rho_ff, Field<VarType::Vector,PhiInterpType::Planes>* E_rho_ff
195 # ifdef EXPLICIT_EM
197 # endif
198  );
199 #else
200 extern "C" void set_rho_pointers(int nrho, int n, double* pot0, double* dpot, Field<VarType::Vector2D,PhiInterpType::None>* E_rho
201 # ifdef SONIC_GK
203 # endif
204  );
205 #endif
206 
207 enum class GridFieldOpts{
208  WithRhoFF = 0,
210 };
211 
212 // Struct for each field that contains different variations of the data: global vs one plane, field-following vs not, maybe host vs device, maybe XGCA vs XGC1
213 template<VarType VT, PhiInterpType PIT>
214 struct GridField{
216 
217  // XGC1
218  Kokkos::View<Field<VT,PIT>*,Kokkos::LayoutRight,HostType> ff_h;
219  Kokkos::View<Field<VT,PIT>**,Kokkos::LayoutRight,HostType> rho_ff_h;
220 
221  GridField(int nphi, int nrho, int nnode, GridFieldOpts opt=GridFieldOpts::WithRhoFF) {
223  ff_h = Kokkos::View<Field<VT,PIT>*,Kokkos::LayoutRight,HostType>("ff_h",nnode);
224  }else{
225  rho_ff_h = Kokkos::View<Field<VT,PIT>**,Kokkos::LayoutRight,HostType>("rho_ff_h",nnode,nrho+1);
226  }
227  }
228 
229  // XGCa
230  Kokkos::View<Field<VT,PIT>**,Kokkos::LayoutRight,HostType> rho_h;
231 
232  // XGCa constructor is distinct since it has 2 arguments vs 3-4 in above XGC1 constructor
233  GridField(int nrho, int nnode)
234  : rho_h("rho_h",nnode, nrho+1) {}
235 };
236 
237 
238 template<class Device, PhiInterpType PIT>
240  private:
243 
244  public:
245 
246  bool turb_efield; // Ideally this would not be here but it is tedious to pass into the gather otherwise
247 
250 
251  // Device Views
252 #ifdef XGC1
253  Kokkos::View<Field<VarType::Vector,PIT>**,Kokkos::LayoutRight,Device> E_phi_ff;
254  Kokkos::View<Field<VarType::Vector,PIT>**,Kokkos::LayoutRight,Device> E_rho_ff;
255  Kokkos::View<Field<VarType::Scalar,PIT>*,Kokkos::LayoutRight,Device> dpot_ff;
256 # ifdef DELTAF_CONV
257  Kokkos::View<Field<VarType::Vector2D,PIT>*,Kokkos::LayoutRight,Device> E00_ff;
258  Kokkos::View<Field<VarType::Scalar,PIT>**,Kokkos::LayoutRight,Device> ddpotdt_phi_ff;
259 # endif
260 # ifdef EXPLICIT_EM
261  //For EM
262  Kokkos::View<Field<VarType::Vector,PIT>**,Kokkos::LayoutRight,Device> dAh_phi_ff;
263  Kokkos::View<Field<VarType::Vector,PIT>**,Kokkos::LayoutRight,Device> dAh_rho_ff;
264  Kokkos::View<Field<VarType::Scalar,PIT>**,Kokkos::LayoutRight,Device> Ah_phi_ff;
265  Kokkos::View<Field<VarType::Scalar,PIT>**,Kokkos::LayoutRight,Device> Ah_rho_ff;
266  Kokkos::View<Field<VarType::Vector,PIT>**,Kokkos::LayoutRight,Device> dAs_phi_ff;
267  Kokkos::View<Field<VarType::Vector,PIT>**,Kokkos::LayoutRight,Device> dAs_rho_ff;
268  Kokkos::View<Field<VarType::Scalar,PIT>**,Kokkos::LayoutRight,Device> As_phi_ff;
269  Kokkos::View<Field<VarType::Scalar,PIT>**,Kokkos::LayoutRight,Device> As_rho_ff;
270  Kokkos::View<Field<VarType::Scalar,PIT>**,Kokkos::LayoutRight,Device> Ah_cv_phi_ff;
271  Kokkos::View<Field<VarType::Scalar,PIT>**,Kokkos::LayoutRight,Device> Epar_em_phi_ff;
272  Kokkos::View<Field<VarType::Scalar,PIT>**,Kokkos::LayoutRight,Device> Epar_em_rho_ff;
273 # endif
274 #else
275  Kokkos::View<Field<VarType::Vector2D,PIT>**,Kokkos::LayoutRight,Device> E_rho;
276  Kokkos::View<double*,Kokkos::LayoutRight,Device> dpot;
277 # ifdef SONIC_GK
278  View<Field<VarType::Vector2D,PIT>**,CLayout,Device> du2_E_rho;
279 # endif
280 #endif
281 
283  : is_initialized(false) {}
284 
285  GridFieldPack(KinType kintype_in, bool turb_efield_in)
286  : is_initialized(true),
287  kintype(kintype_in),
288  turb_efield(turb_efield_in),
289  phi_offset(0),
290  node_offset(0) {}
291 
292  bool can_reuse(KinType kintype_in){
293  if(is_initialized){
294 #ifdef XGC1
295  // For XGC1, can reuse if the existing GridFieldPack is the same KinType
296  return kintype_in==kintype;
297 #else
298  // For XGCa, the drift kinetic GridFieldPack is a subset of the gyrokinetic one,
299  // so is always reused
300  return true;
301 #endif
302  }else{
303  // Cant reuse if uninitialized
304  return false;
305  }
306  }
307 
308  // Device functions:
309 
310  // A correction to the field
311  KOKKOS_INLINE_FUNCTION void field_correction(int i_simd, double basis, double gamma_psi, const SimdVector2D &gradpsi, double* rvec, double* zvec) const {
312  rvec[0]=basis + (1.0-basis) * gamma_psi * gradpsi.r[i_simd];
313  rvec[1]= (1.0-basis) * gamma_psi * gradpsi.z[i_simd];
314  zvec[0]= (1.0-basis) * gamma_psi *(-gradpsi.z[i_simd]);
315  zvec[1]=basis + (1.0-basis) * gamma_psi * gradpsi.r[i_simd];
316  }
317 
318  KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector& vec, int i_simd, const SimdVector& B, double Bmag) const{
319  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];
320  }
321 
322 
323 
324  KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int i_node, double wp, const double (&rvec)[2], const double (&zvec)[2], const FieldWeights<DriftKin, PIT>& wts, LocalFields& fld) const {
325  const int node = i_node - node_offset;
326 #ifdef XGCA
327  E_rho(node,0).gather(fld.E,i_simd, wp, rvec, zvec); // irho = 0
328 #else
329  const int i_phi = (wts.phi.i - phi_offset)%E_phi_ff.extent(0); //grid.nplanes;
330  E_phi_ff(i_phi,node).gather(fld.E,i_simd, wp, wts.phi.w, rvec, zvec);
331 #ifdef EXPLICIT_EM
332  Ah_phi_ff(i_phi,node).gather(fld.Ah,i_simd, wp, wts.phi.w);
333  dAh_phi_ff(i_phi,node).gather(fld.dAh,i_simd, wp, wts.phi.w, rvec, zvec);
334 
335  if(push_controls.em_control_variate) {
336  Ah_cv_phi_ff(i_phi,node).gather(fld.Ah_cv,i_simd, wp, wts.phi.w);
337  }
338 
339  if(push_controls.em_mixed_variable) {
340  dAs_phi_ff(i_phi,node).gather(fld.dAs,i_simd, wp, wts.phi.w, rvec, zvec);
341  As_phi_ff(i_phi,node).gather(fld.As,i_simd, wp, wts.phi.w);
342  }
343  if(push_controls.em_mixed_variable && (push_controls.em_pullback_mode==4)) {
344  Epar_em_phi_ff(i_phi,node).gather(fld.Epar_em,i_simd,wp,wts.phi.w);
345  }
346 #endif
347 #endif
348 
349 #if defined(XGC1) && defined(DELTAF_CONV)
350  E00_ff(node).gather(fld.E00,i_simd, wp, wts.phi.w, rvec, zvec);
351  ddpotdt_phi_ff(i_phi,node).gather(fld.ddpotdt,i_simd, wp, wts.phi.w);
352 #endif
353  }
354 
355  // Gather all fields (ion)
356  KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int node, double wp, const double (&rvec)[2], const double (&zvec)[2], const FieldWeights<GyroKin, PIT>& wts, LocalFields& fld) const {
357 #ifdef NEWGYROMATRIX
358  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
359  int irho = wts.rho[0].i;
360 #else
361  int irho = wts.rho.i;
362 #endif
363 #ifdef XGCA
364  E_rho(node,irho).gather(fld.E,i_simd, wp, wts, rvec, zvec, E_rho(node,irho+1));
365 # ifdef SONIC_GK
366  du2_E_rho(node,0).gather(fld.du2_E,i_simd, wp, rvec, zvec); // irho = 0
367 # endif
368 #else
369  E_rho_ff(node,irho).gather(fld.E,i_simd, wp, wts, rvec, zvec, E_rho_ff(node,irho+1));
370 #ifdef EXPLICIT_EM
371  Ah_rho_ff(node,irho).gather(fld.Ah,i_simd, wp, wts, Ah_rho_ff(node,irho+1));
372  dAh_rho_ff(node,irho).gather(fld.dAh,i_simd, wp, wts, rvec, zvec, dAh_rho_ff(node,irho+1));
373 
374  if(push_controls.em_mixed_variable) {
375  dAs_rho_ff(node,irho).gather(fld.dAs,i_simd, wp, wts, rvec, zvec, dAs_rho_ff(node,irho+1));
376  As_rho_ff(node,irho).gather(fld.As,i_simd, wp, wts, As_rho_ff(node,irho+1));
377  }
378  if(push_controls.em_mixed_variable && (push_controls.em_pullback_mode==4)) {
379  Epar_em_rho_ff(node,irho).gather(fld.Epar_em,i_simd,wp,wts, Epar_em_rho_ff(node,irho+1));
380  }
381 #endif
382 #endif
383  }
384 
385 
397  template<KinType PT>
398  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
399  {
400  // Initialize local fields to 0
401  fld.E.zero();
402 #ifdef DELTAF_CONV
403  fld.E00.zero();
404  fld.ddpotdt.zero();
405 #endif
406 #ifdef EXPLICIT_EM
407  fld.dAh.zero();
408  fld.dAs.zero();
409  fld.Ah.zero();
410  fld.As.zero();
411  fld.Ah_cv.zero();
412  fld.Epar_em.zero();
413 #endif
414 #ifdef SONIC_GK
415  fld.du2_E.zero();
416 #endif
417 
418  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
419  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
420  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
421 
422  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
423 
424  // Get rho and phi weights
425  FieldWeights<PT, PIT> wts(grid, phi, i_simd, rho);
426 
427  double gamma_psi=1.0/sqrt(gradpsi.r[i_simd]*gradpsi.r[i_simd] + gradpsi.z[i_simd]*gradpsi.z[i_simd]);
428 
429  double rvec[2];
430  double zvec[2];
431 
432  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
433  for (int ip = 0; ip<3; ip++){
434  int node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
435  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
436 
437  field_correction(i_simd, grid.basis(node), gamma_psi, gradpsi, rvec, zvec);
438 
439  gather_all_fields(push_controls, i_simd, node, wp, rvec, zvec, wts, fld);
440  }
441 
442 #ifdef NEOCLASSICAL_TEST
443  fld.E.phi[i_simd]=0.0;
444 #else
445 #ifdef XGC1
446  if(turb_efield){
447  // Correct phi component of E-field
448  double Bmag=sqrt(B.r[i_simd]*B.r[i_simd] + B.z[i_simd]*B.z[i_simd] + B.phi[i_simd]*B.phi[i_simd]);
449 
450  // Get phi from para
451  phi_from_para(fld.E, i_simd, B, Bmag);
452 #ifdef EXPLICIT_EM
453  phi_from_para(fld.dAh, i_simd, B, Bmag);
454  phi_from_para(fld.dAs, i_simd, B, Bmag);
455 #endif
456  }else{
457  fld.E.phi[i_simd]=0.0;
458 #ifdef EXPLICIT_EM
459  fld.dAh.phi[i_simd]=0.0;
460  fld.dAs.phi[i_simd]=0.0;
461 #endif
462  }
463 #endif
464 #endif
465  }
466  }
467 
468 
469  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
470  {
471  // Initialize local fields to 0
472 #ifdef EXPLICIT_EM
473  Ah = 0.0;
474 
475  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
476  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
477  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
478  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
479 
480  // Get rho and phi weights
481  FieldWeights<GyroKin, PIT> wts(grid, phi, i_simd, rho);
482 #ifdef NEWGYROMATRIX
483  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
484  int irho = wts.rho[0].i;
485 #else
486  int irho = wts.rho.i;
487 #endif
488 
489  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
490  for (int ip = 0; ip<3; ip++){
491  int node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
492  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
493 
494  Ah_rho_ff(node,irho).gather(Ah,i_simd, wp, wts, Ah_rho_ff(node,irho+1));
495  }
496  }
497 #endif
498  }
499 
500  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
501  {
502  // Initialize local fields to 0
503 #ifdef EXPLICIT_EM
504  Ah = 0.0;
505 
506  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
507  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
508  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
509  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
510 
511  // Get rho and phi weights
512  FieldWeights<DriftKin, PIT> wts(grid, phi, i_simd, rho);
513 
514  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
515  for (int ip = 0; ip<3; ip++){
516  int i_node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
517  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
518 
519  const int i_phi = (wts.phi.i - phi_offset)%grid.nplanes;
520  const int node = i_node - node_offset;
521 
522  Ah_phi_ff(i_phi,node).gather(Ah,i_simd, wp, wts.phi.w);
523  }
524  }
525 #endif
526  }
527 
528  KOKKOS_INLINE_FUNCTION void get_Ah_cv(const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<DriftKin>& rho, Simd<double>& Ah_cv) const
529  {
530  // Initialize local fields to 0
531 #ifdef EXPLICIT_EM
532  Ah_cv = 0.0;
533 
534  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
535  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
536  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
537  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
538 
539  // Get rho and phi weights
540  FieldWeights<DriftKin, PIT> wts(grid, phi, i_simd, rho);
541 
542  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
543  for (int ip = 0; ip<3; ip++){
544  int i_node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
545  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
546 
547  const int i_phi = (wts.phi.i - phi_offset)%grid.nplanes;
548  const int node = i_node - node_offset;
549 
550  Ah_cv_phi_ff(i_phi,node).gather(Ah_cv,i_simd, wp, wts.phi.w);
551  }
552  }
553 #endif
554  }
555 
556  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 {
557 
558  dpot_out.zero();
559 
560  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
561  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
562  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
563 
564  double phi = itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
566  FieldWeights<DriftKin, PIT> wts(grid, phi, i_simd, rho);
567 
568  for(int ip=0; ip<3; ip++){
569  int node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
570  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
571 
572 #ifdef XGC1
573  //#ifndef EXPLICIT_EM
574  // for electron -- rho=0 case, optimized.
575  dpot_ff(node).gather(dpot_out,i_simd, wp, wts.phi.w);
576  //#else
577  // dpot_out[i_simd]+=wp*dpot_n0(node);
578  //#endif
579 #elif defined(XGCA)
580  // Simple gather
581  dpot_out[i_simd]+=wp*dpot(node);
582 #endif
583  }
584  }
585  }
586 };
587 #endif
Simd< double > r
Definition: simd.hpp:145
KOKKOS_INLINE_FUNCTION void gather(SimdVector &vec, int i_simd, double wp, const double(&rvec)[2], const double(&zvec)[2]) const
Definition: grid_field_pack.hpp:154
Definition: simd.hpp:144
GridFieldPack()
Definition: grid_field_pack.hpp:282
void set_rho_pointers(int nrho, int n, double *pot0, double *dpot, Field< VarType::Vector2D, PhiInterpType::None > *E_rho)
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:556
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:469
Kokkos::View< Field< VarType::Vector2D, PIT > **, Kokkos::LayoutRight, Device > E_rho
Definition: grid_field_pack.hpp:275
Definition: grid_field_pack.hpp:107
bool turb_efield
Definition: grid_field_pack.hpp:246
Definition: field_weights.hpp:13
Kokkos::Device< HostExSpace, HostMemSpace > HostType
Definition: space_settings.hpp:56
GridFieldOpts
Definition: grid_field_pack.hpp:207
LinearWeights phi
Definition: field_weights.hpp:35
Definition: push_controls.hpp:8
int phi_offset
Offset for phi_ff field decomposition.
Definition: grid_field_pack.hpp:248
Definition: grid.hpp:10
Kokkos::View< double *, Kokkos::LayoutRight, Device > dpot
Definition: grid_field_pack.hpp:276
int nplanes
Number of planes.
Definition: grid.hpp:84
GridField(int nphi, int nrho, int nnode, GridFieldOpts opt=GridFieldOpts::WithRhoFF)
Definition: grid_field_pack.hpp:221
Definition: grid_field_pack.hpp:239
Kokkos::LayoutRight CLayout
Definition: space_settings.hpp:60
Kokkos::View< Field< VT, PIT > **, Kokkos::LayoutRight, HostType > rho_h
Definition: grid_field_pack.hpp:230
Definition: local_fields.hpp:7
Definition: grid_field_pack.hpp:214
bool can_reuse(KinType kintype_in)
Definition: grid_field_pack.hpp:292
double S[2]
Definition: grid_field_pack.hpp:108
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:500
Definition: field_weights.hpp:33
Definition: grid_field_pack.hpp:27
LinearWeights rho
Definition: field_weights.hpp:36
KOKKOS_INLINE_FUNCTION void get_Ah_cv(const Grid< Device > &grid, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, SimdGyroRadius< DriftKin > &rho, Simd< double > &Ah_cv) const
Definition: grid_field_pack.hpp:528
r coordinate
Definition: globals.hpp:102
KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int i_node, double wp, const double(&rvec)[2], const double(&zvec)[2], const FieldWeights< DriftKin, PIT > &wts, LocalFields &fld) const
Definition: grid_field_pack.hpp:324
KinType kintype
Definition: grid_field_pack.hpp:241
Kokkos::View< int *, Kokkos::LayoutRight, Device > basis
A basis for the guesses?
Definition: grid.hpp:91
Definition: field_weights.hpp:62
KOKKOS_INLINE_FUNCTION void gather(SimdVector &vec, int i_simd, double wp, const double(&wphi)[2], const double(&rvec)[2], const double(&zvec)[2]) const
Definition: grid_field_pack.hpp:97
double S
Definition: grid_field_pack.hpp:180
Simd< double > z
Definition: simd.hpp:141
KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector &vec, int i_simd, const SimdVector &B, double Bmag) const
Definition: grid_field_pack.hpp:318
double w[2]
Definition: linear_weights.hpp:9
KOKKOS_INLINE_FUNCTION void gather(SimdVector &vec, int i_simd, double wp, const FieldWeights< GyroKin, PhiInterpType::Planes > &wts, const double(&rvec)[2], const double(&zvec)[2], const Field< VarType::Vector, PhiInterpType::Planes > &field2) const
Definition: grid_field_pack.hpp:56
KOKKOS_INLINE_FUNCTION void gather(Simd< double > &sca, int i_simd, double wp, const FieldWeights< GyroKin, PhiInterpType::Planes > &wts, const Field< VarType::Scalar, PhiInterpType::Planes > &field2) const
Definition: grid_field_pack.hpp:127
GridFieldPack(KinType kintype_in, bool turb_efield_in)
Definition: grid_field_pack.hpp:285
KOKKOS_INLINE_FUNCTION void gather(Simd< double > &sca, int i_simd, double wp)
Definition: grid_field_pack.hpp:187
Simd< double > phi
Definition: simd.hpp:147
int node_offset
Offset for phi_ff field decomposition.
Definition: grid_field_pack.hpp:249
Simd< double > z
Definition: simd.hpp:146
KinType
Definition: globals.hpp:81
Definition: grid_field_pack.hpp:23
Kokkos::View< Field< VT, PIT > **, Kokkos::LayoutRight, HostType > rho_ff_h
Contains the local plane of gyroaveraged, field-following, host field.
Definition: grid_field_pack.hpp:219
KOKKOS_INLINE_FUNCTION void zero()
Definition: simd.hpp:76
Definition: grid_structs.hpp:7
VarType
Definition: grid_field_pack.hpp:16
Kokkos::View< Field< VT, PIT > *, Kokkos::LayoutRight, HostType > ff_h
Contains the local plane, field-following, host field.
Definition: grid_field_pack.hpp:218
Simd< double > r
Definition: simd.hpp:140
KOKKOS_INLINE_FUNCTION void gather(SimdVector &vec, int i_simd, double wp, const FieldWeights< GyroKin, PhiInterpType::None > &wts, const double(&rvec)[2], const double(&zvec)[2], const Field< VarType::Vector2D, PhiInterpType::None > &field2) const
Definition: grid_field_pack.hpp:168
GridField(int nrho, int nnode)
Definition: grid_field_pack.hpp:233
Definition: simd.hpp:139
KOKKOS_INLINE_FUNCTION void gather(Simd< double > &sca, int i_simd, double wp, const double(&wphi)[2]) const
Definition: grid_field_pack.hpp:116
KOKKOS_INLINE_FUNCTION void zero()
Definition: simd.hpp:162
phi coordinate
Definition: globals.hpp:104
bool is_initialized
Definition: grid_field_pack.hpp:242
Definition: gyro_radius.hpp:96
LinearWeights rho
Definition: field_weights.hpp:64
GridField()
Definition: grid_field_pack.hpp:215
Definition: gyro_radius.hpp:74
KOKKOS_INLINE_FUNCTION void gather(SimdVector &vec, int i_simd, double wp, const double(&wphi)[2], const double(&rvec)[2], const double(&zvec)[2]) const
Definition: grid_field_pack.hpp:38
z coordinate
Definition: globals.hpp:103
double E[2]
Definition: grid_field_pack.hpp:145
KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int node, double wp, const double(&rvec)[2], const double(&zvec)[2], const FieldWeights< GyroKin, PIT > &wts, LocalFields &fld) const
Definition: grid_field_pack.hpp:356
SimdVector E
Definition: local_fields.hpp:8
double V[3][2]
Definition: grid_field_pack.hpp:28
KOKKOS_INLINE_FUNCTION void field_correction(int i_simd, double basis, double gamma_psi, const SimdVector2D &gradpsi, double *rvec, double *zvec) const
Definition: grid_field_pack.hpp:311
Definition: grid_field_pack.hpp:144
Kokkos::View< Vertex *, Kokkos::LayoutRight, Device > nodes
Definition: grid.hpp:92
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:398