XGC1
 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 "gyro_radius.hpp"
12 #include "grid_field.hpp"
13 #include "label.hpp"
14 #include "pack.hpp"
15 
16 // Base struct for GridFieldPack so that non-GPU code doesn't need to be templated just because there is a GridFieldPack argument
18 };
19 
20 template<class Device, PhiInterpType PIT>
22  // Temporary implementation of this to make the PRs more legible
23 #ifdef EXPLICIT_EM
24  using pack_type = Pack<Labeled<EfieldType, Label::E>,
30 
37 #else
38 # ifdef DELTAF_CONV
39  using pack_type = Pack<Labeled<EfieldType, Label::E>,
43 # else
44  using pack_type = Pack<Labeled<EfieldType, Label::E>,
46 # endif
47 #endif
48 
49  private:
52 
53  public:
54 
55  bool turb_efield; // Ideally this would not be here but it is tedious to pass into the gather otherwise
56 
57  int phi_offset;
59 
60 #ifdef XGC1
61  // Views stored in Pack
63 
64  // Const access to pack contents for inside kernels (can modify data but not structure)
65  template<Label FN>
66  KOKKOS_INLINE_FUNCTION const auto& get() const{
67  return pack.template get<FN>();
68  }
69 
70  // Non-const access to pack contents (can modify structure)
71  template<Label FN>
72  inline auto& get(){
73  return pack.template get<FN>();
74  }
75 #endif
76 
77 
78 #ifdef XGC1
80 #else
83 # ifdef SONIC_GK
87 # endif
88 #endif
89  // This will be a ScalarGridField in case of XGCa and a grid field with field-following
90  // interpolation in case of XGC1
92 
94  : is_initialized(false) {}
95 
96  GridFieldPack(KinType kintype_in, bool turb_efield_in)
97  : is_initialized(true),
98  kintype(kintype_in),
99  turb_efield(turb_efield_in),
100  phi_offset(0),
101  node_offset(0) {}
102 
103  bool can_reuse(KinType kintype_in){
104  if(is_initialized){
105 #ifdef XGC1
106  // For XGC1, can reuse if the existing GridFieldPack is the same KinType
107  return kintype_in==kintype;
108 #else
109  // For XGCa, the drift kinetic GridFieldPack is a subset of the gyrokinetic one,
110  // so is always reused
111  return true;
112 #endif
113  }else{
114  // Cant reuse if uninitialized
115  return false;
116  }
117  }
118 
119  // Device functions:
120 
121  KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector& vec, int i_simd, const SimdVector& B, double Bmag) const{
122  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];
123  }
124 
125 
126 
127  KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int i_node, double wp, const FieldCorrection& corr, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, const SimdGyroWeights<DriftKin>& rho_wts, LocalFields& fld) const {
128  const int node = i_node - node_offset;
129 #ifdef XGCA
130  E_rho(node,0).gather(fld.template get<Label::E>(),i_simd, wp, corr); // irho = 0
131  // Loop voltage
132  //loop_voltage(node).gather(fld.loop_voltage,i_simd,wp);
133 #else
134  const int i_phi = (grid_wts.phi_wts.phi.i[i_simd] - phi_offset)%get<Label::E>().f.extent(0); //should be grid.nplanes (?);
135  get<Label::E>()(i_phi,node).gather(fld.template get<Label::E>(),i_simd, wp, grid_wts.phi_wts, corr);
136 #ifdef EXPLICIT_EM
137  get<Label::Ah>()(i_phi,node).gather(fld.template get<Label::Ah>(),i_simd, wp, grid_wts.phi_wts);
138  get<Label::dAh>()(i_phi,node).gather(fld.template get<Label::dAh>(),i_simd, wp, grid_wts.phi_wts, corr);
139 
140  if(push_controls.em_mixed_variable) {
141  get<Label::dAs>()(i_phi,node).gather(fld.template get<Label::dAs>(),i_simd, wp, grid_wts.phi_wts, corr);
142  get<Label::As>()(i_phi,node).gather(fld.template get<Label::As>(),i_simd, wp, grid_wts.phi_wts);
143  }
144  if(push_controls.em_mixed_variable && (push_controls.em_pullback_mode==4)) {
145  get<Label::Epar_em>()(i_phi,node).gather(fld.template get<Label::Epar_em>(),i_simd,wp,grid_wts.phi_wts);
146  }
147 #endif
148  // Loop voltage
149  //loop_voltage(node).gather(fld.loop_voltage,i_simd,wp,grid_wts.phi_wts);
150 #endif
151 
152 #if defined(DELTAF_CONV) && !defined(EXPLICIT_EM)
153  get<Label::E00>()(node).gather(fld.template get<Label::E00>(),i_simd, wp, grid_wts.phi_wts, corr);
154  get<Label::ddpotdt>()(i_phi,node).gather(fld.template get<Label::ddpotdt>(),i_simd, wp, grid_wts.phi_wts);
155 #endif
156  }
157 
158  // Gather all fields (ion)
159  KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int node, double wp, const FieldCorrection& corr, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, const SimdGyroWeights<GyroKin>& rho_wts, LocalFields& fld) const {
160  int irho = rho_wts.irho(i_simd);
161 #ifdef XGCA
162  E_rho(node,irho).gather(fld.template get<Label::E>(),i_simd, wp, rho_wts, corr, E_rho(node,irho+1));
163 # ifdef SONIC_GK
164  dEr_B2_rho(node,0).gather(fld.template get<Label::dEr_B2>(),i_simd, wp, corr); // irho = 0
165  dEz_B2_rho(node,0).gather(fld.template get<Label::dEz_B2>(),i_simd, wp, corr); // irho = 0
166  du2_E_rho(node,0).gather(fld.template get<Label::du2_E>(),i_simd, wp, corr); // irho = 0
167 # endif
168  // Loop voltage
169  //loop_voltage(node).gather(fld.loop_voltage,i_simd,wp);
170 #else
171  get<Label::E_gyro>()(node,irho).gather(fld.template get<Label::E>(),i_simd, wp, rho_wts, grid_wts.phi_wts, corr, get<Label::E_gyro>()(node,irho+1));
172 #ifdef EXPLICIT_EM
173  get<Label::Ah_gyro>()(node,irho).gather(fld.template get<Label::Ah>(),i_simd, wp, rho_wts, grid_wts.phi_wts, get<Label::Ah_gyro>()(node,irho+1));
174  get<Label::dAh_gyro>()(node,irho).gather(fld.template get<Label::dAh>(),i_simd, wp, rho_wts, grid_wts.phi_wts, corr, get<Label::dAh_gyro>()(node,irho+1));
175 
176  if(push_controls.em_mixed_variable) {
177  get<Label::dAs_gyro>()(node,irho).gather(fld.template get<Label::dAs>(),i_simd, wp, rho_wts, grid_wts.phi_wts, corr, get<Label::dAs_gyro>()(node,irho+1));
178  get<Label::As_gyro>()(node,irho).gather(fld.template get<Label::As>(),i_simd, wp, rho_wts, grid_wts.phi_wts, get<Label::As_gyro>()(node,irho+1));
179  }
180  if(push_controls.em_mixed_variable && (push_controls.em_pullback_mode==4)) {
181  get<Label::Epar_em_gyro>()(node,irho).gather(fld.template get<Label::Epar_em>(),i_simd,wp,rho_wts, grid_wts.phi_wts, get<Label::Epar_em_gyro>()(node,irho+1));
182  }
183 #endif
184  // Loop voltage
185  //loop_voltage(node).gather(fld.loop_voltage,i_simd,wp,grid_wts.phi_wts);
186 #endif
187  }
188 
189 
200  template<KinType PT>
201  KOKKOS_INLINE_FUNCTION void fields_at_point(const PushControls &push_controls, const Grid<Device> &grid, const SimdVector &B, const SimdVector2D &gradpsi, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, const SimdGyroWeights<PT>& rho_wts, LocalFields& fld) const
202  {
203  // Initialize local fields to 0
204  fld.template get<Label::E>().zero();
205 #if defined(DELTAF_CONV) && !defined(EXPLICIT_EM)
206  fld.template get<Label::E00>().zero();
207  fld.template get<Label::ddpotdt>().zero();
208 #endif
209 #ifdef EXPLICIT_EM
210  fld.template get<Label::dAh>().zero();
211  fld.template get<Label::dAs>().zero();
212  fld.template get<Label::Ah>().zero();
213  fld.template get<Label::As>().zero();
214  fld.template get<Label::Epar_em>().zero();
215 #endif
216 #ifdef SONIC_GK
217  fld.template get<Label::dEr_B2>()zero();
218  fld.template get<Label::dEz_B2>()zero();
219  fld.template get<Label::du2_E>()zero();
220 #endif
221 
222  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
223  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
224  int itr_work = grid_wts.is_valid(i_simd) ? grid_wts.itr[i_simd] : 1;
225 
226  // Field basis correction
227  FieldCorrection field_correction(gradpsi, i_simd);
228 
229  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
230  for (int ip = 0; ip<3; ip++){
231  int node = grid_wts.is_valid(i_simd) ? grid.get_node_index(itr_work, ip) : 0;
232  double wp = grid_wts.is_valid(i_simd) ? grid_wts.p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
233 
234  field_correction.set(i_simd, grid.uses_rz_basis(node), gradpsi);
235 
236  gather_all_fields(push_controls, i_simd, node, wp, field_correction, grid_wts, rho_wts, fld);
237  }
238 
239 #ifdef NEOCLASSICAL_TEST
240  fld.template get<Label::E>().phi[i_simd]=0.0;
241 #else
242 #ifdef XGC1
243  if(turb_efield){
244  // Correct phi component of E-field
245  double Bmag = B.magnitude(i_simd);
246 
247  // Get phi from para
248  phi_from_para(fld.template get<Label::E>(), i_simd, B, Bmag);
249 #ifdef EXPLICIT_EM
250  phi_from_para(fld.template get<Label::dAh>(), i_simd, B, Bmag);
251  phi_from_para(fld.template get<Label::dAs>(), i_simd, B, Bmag);
252 #endif
253  }else{
254  fld.template get<Label::E>().phi[i_simd]=0.0;
255 #ifdef EXPLICIT_EM
256  fld.template get<Label::dAh>().phi[i_simd]=0.0;
257  fld.template get<Label::dAs>().phi[i_simd]=0.0;
258 #endif
259  }
260 #endif
261 #endif
262  }
263  }
264 
265 
266  KOKKOS_INLINE_FUNCTION void get_Ah(const Grid<Device> &grid, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, const SimdGyroWeights<GyroKin>& rho_wts, Simd<double>& Ah) const
267  {
268  // Initialize local fields to 0
269 #ifdef EXPLICIT_EM
270  Ah = 0.0;
271 
272  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
273  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
274  int itr_work = grid_wts.is_valid(i_simd) ? grid_wts.itr[i_simd] : 1;
275 
276  int irho = rho_wts.irho(i_simd);
277 
278  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
279  for (int ip = 0; ip<3; ip++){
280  int node = grid_wts.is_valid(i_simd) ? grid.get_node_index(itr_work, ip) : 0;
281  double wp = grid_wts.is_valid(i_simd) ? grid_wts.p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
282 
283  get<Label::Ah_gyro>()(node,irho).gather(Ah,i_simd, wp, rho_wts, grid_wts.phi_wts, get<Label::Ah_gyro>()(node,irho+1));
284  }
285  }
286 #endif
287  }
288 
289  KOKKOS_INLINE_FUNCTION void get_Ah(const Grid<Device> &grid, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, const SimdGyroWeights<DriftKin>& rho_wts, Simd<double>& Ah) const
290  {
291  // Initialize local fields to 0
292 #ifdef EXPLICIT_EM
293  Ah = 0.0;
294 
295  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
296  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
297  int itr_work = grid_wts.is_valid(i_simd) ? grid_wts.itr[i_simd] : 1;
298 
299  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
300  for (int ip = 0; ip<3; ip++){
301  int i_node = grid_wts.is_valid(i_simd) ? grid.get_node_index(itr_work, ip) : 0;
302  double wp = grid_wts.is_valid(i_simd) ? grid_wts.p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
303 
304  const int node = i_node - node_offset;
305  constexpr int PHI_ZERO=0;
306  get<Label::Ah>()(PHI_ZERO,node).gather(Ah,i_simd, wp, grid_wts.phi_wts);
307  }
308  }
309 #endif
310  }
311 
312  KOKKOS_INLINE_FUNCTION void get_dpot(const Grid<Device> &grid, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, Simd<double>& dpot_out) const {
313 
314  dpot_out.zero();
315 
316 
317  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
318  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
319  int itr_work = grid_wts.is_valid(i_simd) ? grid_wts.itr[i_simd] : 1;
320 
321  for(int ip=0; ip<3; ip++){
322  int node = grid_wts.is_valid(i_simd) ? grid.get_node_index(itr_work, ip) : 0;
323  double wp = grid_wts.is_valid(i_simd) ? grid_wts.p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
324 
325 #ifdef XGC1
326  //#ifndef EXPLICIT_EM
327  // for electron -- rho=0 case, optimized.
328  dpot_ff(node).gather(dpot_out,i_simd, wp, grid_wts.phi_wts);
329  //#else
330  // dpot_out[i_simd]+=wp*dpot_n0(node);
331  //#endif
332 #elif defined(XGCA)
333  // Simple gather
334  dpot(node).gather(dpot_out,i_simd,wp);
335 #endif
336  }
337  }
338  }
339 };
340 #endif
KOKKOS_INLINE_FUNCTION void set(int i_simd, double basis, const SimdVector2D &gradpsi)
Definition: field.hpp:35
Simd< double > r
Definition: simd.hpp:150
Definition: simd.hpp:149
KOKKOS_INLINE_FUNCTION void get_Ah(const Grid< Device > &grid, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, const SimdGyroWeights< GyroKin > &rho_wts, Simd< double > &Ah) const
Definition: grid_field_pack.hpp:266
GridFieldPack()
Definition: grid_field_pack.hpp:93
GridField< Device, VarType::Scalar, PIT, TorType::OnePlane, KinType::DriftKin > loop_voltage
Definition: grid_field_pack.hpp:91
KOKKOS_INLINE_FUNCTION int irho(int i_simd) const
Definition: gyro_radius.hpp:107
KOKKOS_INLINE_FUNCTION void get_dpot(const Grid< Device > &grid, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, Simd< double > &dpot_out) const
Definition: grid_field_pack.hpp:312
bool turb_efield
Definition: grid_field_pack.hpp:55
KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int i_node, double wp, const FieldCorrection &corr, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, const SimdGyroWeights< DriftKin > &rho_wts, LocalFields &fld) const
Definition: grid_field_pack.hpp:127
Definition: grid_weights.hpp:47
Definition: push_controls.hpp:9
int phi_offset
Offset for phi_ff field decomposition.
Definition: grid_field_pack.hpp:57
Definition: grid.hpp:21
KOKKOS_INLINE_FUNCTION void gather_all_fields(const PushControls &push_controls, int i_simd, int node, double wp, const FieldCorrection &corr, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, const SimdGyroWeights< GyroKin > &rho_wts, LocalFields &fld) const
Definition: grid_field_pack.hpp:159
Definition: grid_field_pack.hpp:21
KOKKOS_INLINE_FUNCTION int uses_rz_basis(const int inode) const
Definition: grid.tpp:248
Definition: local_fields.hpp:33
bool can_reuse(KinType kintype_in)
Definition: grid_field_pack.hpp:103
Definition: gyro_radius.hpp:114
Definition: gyro_radius.hpp:24
Definition: gyro_radius.hpp:84
KinType kintype
Definition: grid_field_pack.hpp:50
KOKKOS_INLINE_FUNCTION int get_node_index(int triangle_index, int tri_vertex_index) const
Definition: grid.tpp:160
Definition: field.hpp:23
KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector &vec, int i_simd, const SimdVector &B, double Bmag) const
Definition: grid_field_pack.hpp:121
Definition: globals.hpp:90
GridFieldPack(KinType kintype_in, bool turb_efield_in)
Definition: grid_field_pack.hpp:96
Definition: label.hpp:44
Simd< double > phi
Definition: simd.hpp:152
int node_offset
Offset for phi_ff field decomposition.
Definition: grid_field_pack.hpp:58
Simd< double > z
Definition: simd.hpp:151
KinType
Definition: globals.hpp:88
KOKKOS_INLINE_FUNCTION void zero()
Definition: simd.hpp:76
KOKKOS_INLINE_FUNCTION void get_Ah(const Grid< Device > &grid, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, const SimdGyroWeights< DriftKin > &rho_wts, Simd< double > &Ah) const
Definition: grid_field_pack.hpp:289
pack_type pack
Definition: grid_field_pack.hpp:62
GridField< Device, VarType::Scalar, PIT, TorType::OnePlane, KinType::DriftKin > dpot_ff
Definition: grid_field_pack.hpp:79
Definition: simd.hpp:139
KOKKOS_INLINE_FUNCTION void fields_at_point(const PushControls &push_controls, const Grid< Device > &grid, const SimdVector &B, const SimdVector2D &gradpsi, const SimdGridWeights< Order::One, PIT_GLOBAL > &grid_wts, const SimdGyroWeights< PT > &rho_wts, LocalFields &fld) const
Definition: grid_field_pack.hpp:201
bool is_initialized
Definition: grid_field_pack.hpp:51
Pack< Labeled< EfieldType, Label::E >, Labeled< EfieldGyroType, Label::E_gyro >> pack_type
Definition: grid_field_pack.hpp:45
KOKKOS_INLINE_FUNCTION double magnitude(const int i_simd) const
Definition: simd.hpp:176
Definition: grid_field_pack.hpp:17