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