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  // This will be a ScalarGridField in case of XGCa and a grid field with field-following
67  // interpolation in case of XGC1
69 
71  : is_initialized(false) {}
72 
73  GridFieldPack(KinType kintype_in, bool turb_efield_in, int nrho_in, double inv_drho_in)
74  : is_initialized(true),
75  kintype(kintype_in),
76  turb_efield(turb_efield_in),
77  nrho(nrho_in),
78  inv_drho(inv_drho_in),
79  phi_offset(0),
80  node_offset(0) {}
81 
82  bool can_reuse(KinType kintype_in){
83  if(is_initialized){
84 #ifdef XGC1
85  // For XGC1, can reuse if the existing GridFieldPack is the same KinType
86  return kintype_in==kintype;
87 #else
88  // For XGCa, the drift kinetic GridFieldPack is a subset of the gyrokinetic one,
89  // so is always reused
90  return true;
91 #endif
92  }else{
93  // Cant reuse if uninitialized
94  return false;
95  }
96  }
97 
98  // Device functions:
99 
100  KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector& vec, int i_simd, const SimdVector& B, double Bmag) const{
101  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];
102  }
103 
104 
105 
106  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 {
107  const int node = i_node - node_offset;
108 #ifdef XGCA
109  E_rho(node,0).gather(fld.E,i_simd, wp, corr); // irho = 0
110  // Loop voltage
111  //loop_voltage(node).gather(fld.loop_voltage,i_simd,wp);
112 #else
113  const int i_phi = (wts.phi.i - phi_offset)%E_phi_ff.f.extent(0); //should be grid.nplanes (?);
114  E_phi_ff(i_phi,node).gather(fld.E,i_simd, wp, wts.phi.w, corr);
115 #ifdef EXPLICIT_EM
116  Ah_phi_ff(i_phi,node).gather(fld.Ah,i_simd, wp, wts.phi.w);
117  dAh_phi_ff(i_phi,node).gather(fld.dAh,i_simd, wp, wts.phi.w, corr);
118 
119  if(push_controls.em_mixed_variable) {
120  dAs_phi_ff(i_phi,node).gather(fld.dAs,i_simd, wp, wts.phi.w, corr);
121  As_phi_ff(i_phi,node).gather(fld.As,i_simd, wp, wts.phi.w);
122  }
123  if(push_controls.em_mixed_variable && (push_controls.em_pullback_mode==4)) {
124  Epar_em_phi_ff(i_phi,node).gather(fld.Epar_em,i_simd,wp,wts.phi.w);
125  }
126 #endif
127  // Loop voltage
128  //loop_voltage(node).gather(fld.loop_voltage,i_simd,wp,wts.phi.w);
129 #endif
130 
131 #if defined(XGC1) && defined(DELTAF_CONV)
132  E00_ff(node).gather(fld.E00,i_simd, wp, wts.phi.w, corr);
133  ddpotdt_phi_ff(i_phi,node).gather(fld.ddpotdt,i_simd, wp, wts.phi.w);
134 #endif
135  }
136 
137  // Gather all fields (ion)
138  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 {
139 #ifdef NEWGYROMATRIX
140  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
141  int irho = wts.rho[0].i;
142 #else
143  int irho = wts.rho.i;
144 #endif
145 #ifdef XGCA
146  E_rho(node,irho).gather(fld.E,i_simd, wp, wts, corr, E_rho(node,irho+1));
147 # ifdef SONIC_GK
148  dEr_B2_rho(node,0).gather(fld.dEr_B2,i_simd, wp, corr); // irho = 0
149  dEz_B2_rho(node,0).gather(fld.dEz_B2,i_simd, wp, corr); // irho = 0
150  du2_E_rho(node,0).gather(fld.du2_E,i_simd, wp, corr); // irho = 0
151 # endif
152  // Loop voltage
153  //loop_voltage(node).gather(fld.loop_voltage,i_simd,wp);
154 #else
155  E_rho_ff(node,irho).gather(fld.E,i_simd, wp, wts, corr, E_rho_ff(node,irho+1));
156 #ifdef EXPLICIT_EM
157  Ah_rho_ff(node,irho).gather(fld.Ah,i_simd, wp, wts, Ah_rho_ff(node,irho+1));
158  dAh_rho_ff(node,irho).gather(fld.dAh,i_simd, wp, wts, corr, dAh_rho_ff(node,irho+1));
159 
160  if(push_controls.em_mixed_variable) {
161  dAs_rho_ff(node,irho).gather(fld.dAs,i_simd, wp, wts, corr, dAs_rho_ff(node,irho+1));
162  As_rho_ff(node,irho).gather(fld.As,i_simd, wp, wts, As_rho_ff(node,irho+1));
163  }
164  if(push_controls.em_mixed_variable && (push_controls.em_pullback_mode==4)) {
165  Epar_em_rho_ff(node,irho).gather(fld.Epar_em,i_simd,wp,wts, Epar_em_rho_ff(node,irho+1));
166  }
167 #endif
168  // Loop voltage
169  //loop_voltage(node).gather(fld.loop_voltage,i_simd,wp,wts.phi.w);
170 #endif
171  }
172 
173 
185  template<KinType PT>
186  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 SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, SimdGyroRadius<PT>& rho, LocalFields& fld) const
187  {
188  // Initialize local fields to 0
189  fld.E.zero();
190 #ifdef DELTAF_CONV
191  fld.E00.zero();
192  fld.ddpotdt.zero();
193 #endif
194 #ifdef EXPLICIT_EM
195  fld.dAh.zero();
196  fld.dAs.zero();
197  fld.Ah.zero();
198  fld.As.zero();
199  fld.Epar_em.zero();
200 #endif
201 #ifdef SONIC_GK
202  fld.dEr_B2.zero();
203  fld.dEz_B2.zero();
204  fld.du2_E.zero();
205 #endif
206 
207  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
208  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
209  int itr_work = grid_wts.is_valid(i_simd) ? grid_wts.itr[i_simd] : 1;
210 
211  double phi= grid_wts.is_valid(i_simd) ? fld_phi[i_simd] : 0.0;
212 
213  // Get rho and phi weights
214  FieldWeights<PT, PIT> wts(grid, phi, i_simd, rho, inv_drho, nrho);
215 
216  // Field basis correction
217  FieldCorrection field_correction(gradpsi, i_simd);
218 
219  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
220  for (int ip = 0; ip<3; ip++){
221  int node = grid_wts.is_valid(i_simd) ? grid.get_node_index(itr_work, ip) : 0;
222  double wp = grid_wts.is_valid(i_simd) ? grid_wts.p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
223 
224  field_correction.set(i_simd, grid.uses_rz_basis(node), gradpsi);
225 
226  gather_all_fields(push_controls, i_simd, node, wp, field_correction, wts, fld);
227  }
228 
229 #ifdef NEOCLASSICAL_TEST
230  fld.E.phi[i_simd]=0.0;
231 #else
232 #ifdef XGC1
233  if(turb_efield){
234  // Correct phi component of E-field
235  double Bmag = B.magnitude(i_simd);
236 
237  // Get phi from para
238  phi_from_para(fld.E, i_simd, B, Bmag);
239 #ifdef EXPLICIT_EM
240  phi_from_para(fld.dAh, i_simd, B, Bmag);
241  phi_from_para(fld.dAs, i_simd, B, Bmag);
242 #endif
243  }else{
244  fld.E.phi[i_simd]=0.0;
245 #ifdef EXPLICIT_EM
246  fld.dAh.phi[i_simd]=0.0;
247  fld.dAs.phi[i_simd]=0.0;
248 #endif
249  }
250 #endif
251 #endif
252  }
253  }
254 
255 
256  KOKKOS_INLINE_FUNCTION void get_Ah(const Grid<Device> &grid, const Simd<double>& fld_phi, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, SimdGyroRadius<GyroKin>& rho, Simd<double>& Ah) const
257  {
258  // Initialize local fields to 0
259 #ifdef EXPLICIT_EM
260  Ah = 0.0;
261 
262  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
263  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
264  int itr_work = grid_wts.is_valid(i_simd) ? grid_wts.itr[i_simd] : 1;
265  double phi= grid_wts.is_valid(i_simd) ? fld_phi[i_simd] : 0.0;
266 
267  // Get rho and phi weights
268  FieldWeights<GyroKin, PIT> wts(grid, phi, i_simd, rho, inv_drho, nrho);
269 #ifdef NEWGYROMATRIX
270  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
271  int irho = wts.rho[0].i;
272 #else
273  int irho = wts.rho.i;
274 #endif
275 
276  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
277  for (int ip = 0; ip<3; ip++){
278  int node = grid_wts.is_valid(i_simd) ? grid.get_node_index(itr_work, ip) : 0;
279  double wp = grid_wts.is_valid(i_simd) ? grid_wts.p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
280 
281  Ah_rho_ff(node,irho).gather(Ah,i_simd, wp, wts, Ah_rho_ff(node,irho+1));
282  }
283  }
284 #endif
285  }
286 
287  KOKKOS_INLINE_FUNCTION void get_Ah(const Grid<Device> &grid, const Simd<double>& fld_phi, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, SimdGyroRadius<DriftKin>& rho, Simd<double>& Ah) const
288  {
289  // Initialize local fields to 0
290 #ifdef EXPLICIT_EM
291  Ah = 0.0;
292 
293  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
294  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
295  int itr_work = grid_wts.is_valid(i_simd) ? grid_wts.itr[i_simd] : 1;
296  double phi= grid_wts.is_valid(i_simd) ? fld_phi[i_simd] : 0.0;
297 
298  // Get rho and phi weights
299  FieldWeights<DriftKin, PIT> wts(grid, phi, i_simd, rho, inv_drho, nrho);
300 
301  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
302  for (int ip = 0; ip<3; ip++){
303  int i_node = grid_wts.is_valid(i_simd) ? grid.get_node_index(itr_work, ip) : 0;
304  double wp = grid_wts.is_valid(i_simd) ? grid_wts.p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
305 
306  const int i_phi = (wts.phi.i - phi_offset)%grid.nplanes;
307  const int node = i_node - node_offset;
308 
309  Ah_phi_ff(i_phi,node).gather(Ah,i_simd, wp, wts.phi.w);
310  }
311  }
312 #endif
313  }
314 
315  KOKKOS_INLINE_FUNCTION void get_dpot(const Grid<Device> &grid, const Simd<double>& fld_phi, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, Simd<double>& dpot_out) const {
316 
317  dpot_out.zero();
318 
319  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
320  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
321  int itr_work = grid_wts.is_valid(i_simd) ? grid_wts.itr[i_simd] : 1;
322 
323  double phi = grid_wts.is_valid(i_simd) ? fld_phi[i_simd] : 0.0;
325  FieldWeights<DriftKin, PIT> wts(grid, phi, i_simd, rho, inv_drho, nrho);
326 
327  for(int ip=0; ip<3; ip++){
328  int node = grid_wts.is_valid(i_simd) ? grid.get_node_index(itr_work, ip) : 0;
329  double wp = grid_wts.is_valid(i_simd) ? grid_wts.p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
330 
331 #ifdef XGC1
332  //#ifndef EXPLICIT_EM
333  // for electron -- rho=0 case, optimized.
334  dpot_ff(node).gather(dpot_out,i_simd, wp, wts.phi.w);
335  //#else
336  // dpot_out[i_simd]+=wp*dpot_n0(node);
337  //#endif
338 #elif defined(XGCA)
339  // Simple gather
340  dpot(node).gather(dpot_out,i_simd,wp);
341 #endif
342  }
343  }
344  }
345 };
346 #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:70
GridField< Device, VarType::Scalar, PIT, TorType::OnePlane, KinType::DriftKin > loop_voltage
Definition: grid_field_pack.hpp:68
Definition: gyro_radius.hpp:24
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
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 SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, SimdGyroRadius< PT > &rho, LocalFields &fld) const
Definition: grid_field_pack.hpp:186
Definition: globals.hpp:89
Definition: grid_weights.hpp:51
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:67
int nplanes
Number of planes.
Definition: grid.hpp:269
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:106
KOKKOS_INLINE_FUNCTION int uses_rz_basis(const int inode) const
Definition: grid.tpp:1008
Definition: local_fields.hpp:7
Definition: grid_field.hpp:25
bool can_reuse(KinType kintype_in)
Definition: grid_field_pack.hpp:82
GridFieldPack(KinType kintype_in, bool turb_efield_in, int nrho_in, double inv_drho_in)
Definition: grid_field_pack.hpp:73
KOKKOS_INLINE_FUNCTION void get_Ah(const Grid< Device > &grid, const Simd< double > &fld_phi, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, SimdGyroRadius< DriftKin > &rho, Simd< double > &Ah) const
Definition: grid_field_pack.hpp:287
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:885
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:100
Definition: globals.hpp:90
KOKKOS_INLINE_FUNCTION void get_dpot(const Grid< Device > &grid, const Simd< double > &fld_phi, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, Simd< double > &dpot_out) const
Definition: grid_field_pack.hpp:315
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:88
KOKKOS_INLINE_FUNCTION void zero()
Definition: simd.hpp:76
GridField< Device, VarType::Scalar, PIT, TorType::OnePlane, KinType::DriftKin > dpot
Definition: grid_field_pack.hpp:59
KOKKOS_INLINE_FUNCTION void get_Ah(const Grid< Device > &grid, const Simd< double > &fld_phi, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, SimdGyroRadius< GyroKin > &rho, Simd< double > &Ah) const
Definition: grid_field_pack.hpp:256
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:138
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