XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
field_gather.hpp
Go to the documentation of this file.
1 #ifndef FIELD_GATHER_HPP
2 #define FIELD_GATHER_HPP
3 
4 #include "globals.hpp"
5 #include "sml.hpp"
6 #include "grid.hpp"
7 #include "field_weights.hpp"
8 #include "gyro_radius.hpp"
9 #include "electric_field.hpp"
10 
11 
12 // A correction to the field
13 KOKKOS_INLINE_FUNCTION void field_correction(int i_simd, double basis, double gamma_psi, const SimdVector2D &gradpsi, double* rvec, double* zvec){
14  rvec[0]=basis + (1.0-basis) * gamma_psi * gradpsi.r[i_simd];
15  rvec[1]= (1.0-basis) * gamma_psi * gradpsi.z[i_simd];
16  zvec[0]= (1.0-basis) * gamma_psi *(-gradpsi.z[i_simd]);
17  zvec[1]=basis + (1.0-basis) * gamma_psi * gradpsi.r[i_simd];
18 }
19 
20 // Struct that contains functions for the field gather that need to be templated (drift vs gyrokinetic)
21 // Generic declaration:
22 template<KinType PT>
23 struct FieldGather {
24 
25  // Gather operations for each field
26 
27  // phi_ff
28  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 VectorFieldPlanes& field) const;
29 
30  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 Vector2DFieldPlanes& field) const;
31 
32  KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp, const double (&wphi)[2], const ScalarFieldPlanes& field) const;
33 
34  // rho_ff
35  KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const FieldWeights<PT>& wts, const double (&rvec)[2], const double (&zvec)[2], const InvVectorFieldPlanes& field, const InvVectorFieldPlanes& field2) const;
36 
37  KOKKOS_INLINE_FUNCTION void gather(Simd<double>& sca, int i_simd, double wp, const FieldWeights<PT>& wts, const ScalarFieldPlanes& field, const ScalarFieldPlanes& field2) const;
38 
39  // axisymmetric
40  KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const double (&rvec)[2], const double (&zvec)[2], const FieldXGCa& field) const;
41 
42  KOKKOS_INLINE_FUNCTION void gather(SimdVector& vec, int i_simd, double wp, const FieldWeights<PT>& wts, const double (&rvec)[2], const double (&zvec)[2], const FieldXGCa& field, const FieldXGCa& field2) const;
43 
44  KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector& vec, int i_simd, const SimdVector& B, double Bmag) const;
45 
46 
47 
48  // Get Ah and Ah_cv
49  template<class Device>
50  KOKKOS_INLINE_FUNCTION void get_Ah(const ElectricField<Device> fields, const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<PT>& rho, Simd<double>& Ah) const;
51 
52  template<class Device>
53  KOKKOS_INLINE_FUNCTION void get_Ah_cv(const ElectricField<Device> fields, const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<PT>& rho, Simd<double>& Ah_cv) const;
54 
55  // Get dpot
56  template<class Device>
57  KOKKOS_INLINE_FUNCTION void get_dpot(const ElectricField<Device>& fields, const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, Simd<double>& dpot_out) const;
58 
59  // A generic version of the function which is specialized below
60  template<class Device, typename T>
61  KOKKOS_INLINE_FUNCTION void gather_all_fields(const ElectricField<Device>& fields, const Simulation<Device> &sml, int i_simd, int node, double wp,
62  const double (&rvec)[2], const double (&zvec)[2], T& wts, LocalFields& fld) const {}
63 
64  // The generic function (formally gk/gk_elec)
65  template<class Device>
66  KOKKOS_INLINE_FUNCTION void fields_at_point(const ElectricField<Device> fields, const Simulation<Device> &sml, const Grid<Device> &grid,
67  const SimdVector &B, const SimdVector2D &gradpsi, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<PT>& rho, LocalFields& fld) const;
68 };
69 
70 #ifdef XGC1
71 
80 template<KinType PT>
81 KOKKOS_INLINE_FUNCTION void FieldGather<PT>::gather(SimdVector& vec, int i_simd, double wp, const double (&wphi)[2], const double (&rvec)[2], const double (&zvec)[2], const VectorFieldPlanes& field) const{
82  vec.r[i_simd] += wp*(wphi[0]*( rvec[0]*field.V[PIR][0] + zvec[0]*field.V[PIZ][0] )
83  + wphi[1]*( rvec[0]*field.V[PIR][1] + zvec[0]*field.V[PIZ][1] ) );
84  vec.z[i_simd] += wp*(wphi[0]*( rvec[1]*field.V[PIR][0] + zvec[1]*field.V[PIZ][0] )
85  + wphi[1]*( rvec[1]*field.V[PIR][1] + zvec[1]*field.V[PIZ][1] ) );
86  vec.phi[i_simd] += wp*(wphi[0]* field.V[PIP][0]
87  + wphi[1]* field.V[PIP][1] );
88 }
89 
99 template<KinType PT>
100 KOKKOS_INLINE_FUNCTION void FieldGather<PT>::gather(SimdVector& vec, int i_simd, double wp, const double (&wphi)[2], const double (&rvec)[2], const double (&zvec)[2], const Vector2DFieldPlanes& field) const{
101  vec.r[i_simd] += wp*(wphi[0]*( rvec[0]*field.V[PIR][0] + zvec[0]*field.V[PIZ][0] )
102  + wphi[1]*( rvec[0]*field.V[PIR][1] + zvec[0]*field.V[PIZ][1] ) );
103  vec.z[i_simd] += wp*(wphi[0]*( rvec[1]*field.V[PIR][0] + zvec[1]*field.V[PIZ][0] )
104  + wphi[1]*( rvec[1]*field.V[PIR][1] + zvec[1]*field.V[PIZ][1] ) );
105 }
106 
114 template<KinType PT>
115 KOKKOS_INLINE_FUNCTION void FieldGather<PT>::gather(Simd<double>& sca, int i_simd, double wp, const double (&wphi)[2], const ScalarFieldPlanes& field) const{
116  sca[i_simd] += wp* ( wphi[0]*field.S[0] + wphi[1]*field.S[1] );
117 }
118 
128 template<>
129 KOKKOS_INLINE_FUNCTION void FieldGather<GyroKin>::gather(SimdVector& vec, int i_simd, double wp, const FieldWeights<GyroKin>& wts, const double (&rvec)[2], const double (&zvec)[2], const InvVectorFieldPlanes& field, const InvVectorFieldPlanes& field2) const{
130 #ifdef NEWGYROMATRIX
131  vec.r[i_simd] += wp*(wts.phi.w[0]*( rvec[0]*(wts.rho[0].w[0]*field.V[0][PIR] + wts.rho[0].w[1]*field2.V[0][PIR])
132  + zvec[0]*(wts.rho[0].w[0]*field.V[0][PIZ] + wts.rho[0].w[1]*field2.V[0][PIZ]) )
133  + wts.phi.w[1]*( rvec[0]*(wts.rho[1].w[0]*field.V[1][PIR] + wts.rho[1].w[1]*field2.V[1][PIR])
134  + zvec[0]*(wts.rho[1].w[0]*field.V[1][PIZ] + wts.rho[1].w[1]*field2.V[1][PIZ]) ) );
135  vec.z[i_simd] += wp*(wts.phi.w[0]*( rvec[1]*(wts.rho[0].w[0]*field.V[0][PIR] + wts.rho[0].w[1]*field2.V[0][PIR])
136  + zvec[1]*(wts.rho[0].w[0]*field.V[0][PIZ] + wts.rho[0].w[1]*field2.V[0][PIZ]) )
137  + wts.phi.w[1]*( rvec[1]*(wts.rho[1].w[0]*field.V[1][PIR] + wts.rho[1].w[1]*field2.V[1][PIR])
138  + zvec[1]*(wts.rho[1].w[0]*field.V[1][PIZ] + wts.rho[1].w[1]*field2.V[1][PIZ]) ) );
139  vec.phi[i_simd] += wp*(wts.phi.w[0]* (wts.rho[0].w[0]*field.V[0][PIP] + wts.rho[0].w[1]*field2.V[0][PIP])
140  + wts.phi.w[1]* (wts.rho[1].w[0]*field.V[1][PIP] + wts.rho[1].w[1]*field2.V[1][PIP]) );
141 #else
142  vec.r[i_simd] += wp*(wts.phi.w[0]*( rvec[0]*(wts.rho.w[0]*field.V[0][PIR] + wts.rho.w[1]*field2.V[0][PIR])
143  + zvec[0]*(wts.rho.w[0]*field.V[0][PIZ] + wts.rho.w[1]*field2.V[0][PIZ]) )
144  + wts.phi.w[1]*( rvec[0]*(wts.rho.w[0]*field.V[1][PIR] + wts.rho.w[1]*field2.V[1][PIR])
145  + zvec[0]*(wts.rho.w[0]*field.V[1][PIZ] + wts.rho.w[1]*field2.V[1][PIZ]) ) );
146  vec.z[i_simd] += wp*(wts.phi.w[0]*( rvec[1]*(wts.rho.w[0]*field.V[0][PIR] + wts.rho.w[1]*field2.V[0][PIR])
147  + zvec[1]*(wts.rho.w[0]*field.V[0][PIZ] + wts.rho.w[1]*field2.V[0][PIZ]) )
148  + wts.phi.w[1]*( rvec[1]*(wts.rho.w[0]*field.V[1][PIR] + wts.rho.w[1]*field2.V[1][PIR])
149  + zvec[1]*(wts.rho.w[0]*field.V[1][PIZ] + wts.rho.w[1]*field2.V[1][PIZ]) ) );
150  vec.phi[i_simd] += wp*(wts.phi.w[0]* (wts.rho.w[0]*field.V[0][PIP] + wts.rho.w[1]*field2.V[0][PIP])
151  + wts.phi.w[1]* (wts.rho.w[0]*field.V[1][PIP] + wts.rho.w[1]*field2.V[1][PIP]) );
152 #endif
153 }
154 
162 template<>
163 KOKKOS_INLINE_FUNCTION void FieldGather<GyroKin>::gather(Simd<double>& sca, int i_simd, double wp, const FieldWeights<GyroKin>& wts, const ScalarFieldPlanes& field, const ScalarFieldPlanes& field2) const{
164 #ifdef NEWGYROMATRIX
165  sca[i_simd] += wp* ( wts.phi.w[0]*wts.rho[0].w[0]*field.S[0]
166  + wts.phi.w[1]*wts.rho[1].w[0]*field.S[1]
167  + wts.phi.w[0]*wts.rho[0].w[1]*field2.S[0]
168  + wts.phi.w[1]*wts.rho[1].w[1]*field2.S[1]);
169 #else
170  sca[i_simd] += wp* ( wts.phi.w[0]*wts.rho.w[0]*field.S[0]
171  + wts.phi.w[1]*wts.rho.w[0]*field.S[1]
172  + wts.phi.w[0]*wts.rho.w[1]*field2.S[0]
173  + wts.phi.w[1]*wts.rho.w[1]*field2.S[1]);
174 #endif
175 }
176 
177 #endif
178 
187 template<KinType PT>
188 KOKKOS_INLINE_FUNCTION void FieldGather<PT>::gather(SimdVector& vec, int i_simd, double wp, const double (&rvec)[2], const double (&zvec)[2], const FieldXGCa& field) const{
189  vec.r[i_simd] += wp*( rvec[0]*field.E[PIR] + zvec[0]*field.E[PIZ] );
190  vec.z[i_simd] += wp*( rvec[1]*field.E[PIR] + zvec[1]*field.E[PIZ] );
191 }
192 
193 template<KinType PT>
194 KOKKOS_INLINE_FUNCTION void FieldGather<PT>::gather(SimdVector& vec, int i_simd, double wp, const FieldWeights<PT>& wts, const double (&rvec)[2], const double (&zvec)[2], const FieldXGCa& field, const FieldXGCa& field2) const{
195  vec.r[i_simd] += wp*( rvec[0]*(wts.rho.w[0]*field.E[PIR] + wts.rho.w[1]*field2.E[PIR])
196  +zvec[0]*(wts.rho.w[0]*field.E[PIZ] + wts.rho.w[1]*field2.E[PIZ]) );
197  vec.z[i_simd] += wp*( rvec[1]*(wts.rho.w[0]*field.E[PIR] + wts.rho.w[1]*field2.E[PIR])
198  +zvec[1]*(wts.rho.w[0]*field.E[PIZ] + wts.rho.w[1]*field2.E[PIZ]) );
199 }
200 
201 
202 template<KinType PT>
203 KOKKOS_INLINE_FUNCTION void FieldGather<PT>::phi_from_para(SimdVector& vec, int i_simd, const SimdVector& B, double Bmag) const{
204  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];
205 }
206 
207 
208 
209 template<>
210 template<class Device, typename T>
211 KOKKOS_INLINE_FUNCTION void FieldGather<DriftKin>::gather_all_fields(const ElectricField<Device>& fields, const Simulation<Device> &sml, int i_simd, int node, double wp, const double (&rvec)[2], const double (&zvec)[2], T& wts, LocalFields& fld) const {
212 
213 #ifdef XGCA
214  gather(fld.E,i_simd, wp, rvec, zvec, fields.E_rho(node,0)); // irho = 0
215 #else
216  gather(fld.E,i_simd, wp, wts.phi.w, rvec, zvec, fields.E_phi_ff(wts.phi.i,node));
217 #ifdef EXPLICIT_EM
218  gather(fld.Ah,i_simd, wp, wts.phi.w, fields.Ah_phi_ff(wts.phi.i,node));
219  gather(fld.dAh,i_simd, wp, wts.phi.w, rvec, zvec, fields.dAh_phi_ff(wts.phi.i,node));
220 
221  if(sml.em_control_variate) {
222  gather(fld.Ah_cv,i_simd, wp, wts.phi.w, fields.Ah_cv_phi_ff(wts.phi.i,node));
223  }
224 
225  if(sml.em_mixed_variable) {
226  gather(fld.dAs,i_simd, wp, wts.phi.w, rvec, zvec, fields.dAs_phi_ff(wts.phi.i,node));
227  gather(fld.As,i_simd, wp, wts.phi.w, fields.As_phi_ff(wts.phi.i,node));
228  }
229  if(sml.em_mixed_variable && (sml.em_pullback_mode==4)) {
230  gather(fld.Epar_em,i_simd,wp,wts.phi.w,fields.Epar_em_phi_ff(wts.phi.i,node));
231  }
232 #endif
233 #endif
234 
235 #if defined(XGC1) && defined(DELTAF_CONV)
236  gather(fld.E00,i_simd, wp, wts.phi.w, rvec, zvec, fields.E00_ff(node));
237  gather(fld.ddpotdt,i_simd, wp, wts.phi.w, fields.ddpotdt_phi(wts.phi.i,node));
238 #endif
239 }
240 
241 // Gather all fields (ion)
242 template<>
243 template<class Device, typename T>
244 KOKKOS_INLINE_FUNCTION void FieldGather<GyroKin>::gather_all_fields(const ElectricField<Device>& fields, const Simulation<Device> &sml, int i_simd, int node, double wp, const double (&rvec)[2], const double (&zvec)[2], T& wts, LocalFields& fld) const {
245 #ifdef NEWGYROMATRIX
246  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
247  int irho = wts.rho[0].i;
248 #else
249  int irho = wts.rho.i;
250 #endif
251 #ifdef XGCA
252  gather(fld.E,i_simd, wp, wts, rvec, zvec, fields.E_rho(node,irho), fields.E_rho(node,irho+1));
253 #else
254  gather(fld.E,i_simd, wp, wts, rvec, zvec, fields.E_rho_ff(node,irho),fields.E_rho_ff(node,irho+1));
255 #ifdef EXPLICIT_EM
256  gather(fld.Ah,i_simd, wp, wts, fields.Ah_rho_ff(node,irho),fields.Ah_rho_ff(node,irho+1));
257  gather(fld.dAh,i_simd, wp, wts, rvec, zvec, fields.dAh_rho_ff(node,irho),fields.dAh_rho_ff(node,irho+1));
258 
259  if(sml.em_mixed_variable) {
260  gather(fld.dAs,i_simd, wp, wts, rvec, zvec, fields.dAs_rho_ff(node,irho),fields.dAs_rho_ff(node,irho+1));
261  gather(fld.As,i_simd, wp, wts, fields.As_rho_ff(node,irho), fields.As_rho_ff(node,irho+1));
262  }
263  if(sml.em_mixed_variable && (sml.em_pullback_mode==4)) {
264  gather(fld.Epar_em,i_simd,wp,wts, fields.Epar_em_rho_ff(node,irho), fields.Epar_em_rho_ff(node,irho+1));
265  }
266 #endif
267 #endif
268 }
269 
270 
283 template<KinType PT>
284 template<class Device>
285 KOKKOS_INLINE_FUNCTION void FieldGather<PT>::fields_at_point(const ElectricField<Device> fields, const Simulation<Device> &sml, 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
286 {
287  // Initialize local fields to 0
288  fld.E.zero();
289 #ifdef DELTAF_CONV
290  fld.E00.zero();
291  fld.ddpotdt.zero();
292 #endif
293 #ifdef EXPLICIT_EM
294  fld.dAh.zero();
295  fld.dAs.zero();
296  fld.Ah.zero();
297  fld.As.zero();
298  fld.Ah_cv.zero();
299  fld.Epar_em.zero();
300 #endif
301 
302  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
303  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
304  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
305 
306  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
307 
308  // Get rho and phi weights
309  FieldWeights<PT> wts(grid, phi, i_simd, rho);
310 
311  double gamma_psi=1.0/sqrt(gradpsi.r[i_simd]*gradpsi.r[i_simd] + gradpsi.z[i_simd]*gradpsi.z[i_simd]);
312 
313  double rvec[2];
314  double zvec[2];
315 
316  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
317  for (int ip = 0; ip<3; ip++){
318  int node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
319  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
320 
321  field_correction(i_simd, grid.basis(node), gamma_psi, gradpsi, rvec, zvec);
322 
323  gather_all_fields(fields, sml, i_simd, node, wp, rvec, zvec, wts, fld);
324  }
325 
326 #ifdef NEOCLASSICAL_TEST
327  fld.E.phi[i_simd]=0.0;
328 #else
329 #ifdef XGC1
330  if(fields.turb_efield){
331  // Correct phi component of E-field
332  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]);
333 
334  // Get phi from para
335  phi_from_para(fld.E, i_simd, B, Bmag);
336 #ifdef EXPLICIT_EM
337  phi_from_para(fld.dAh, i_simd, B, Bmag);
338  phi_from_para(fld.dAs, i_simd, B, Bmag);
339 #endif
340  }else{
341  fld.E.phi[i_simd]=0.0;
342 #ifdef EXPLICIT_EM
343  fld.dAh.phi[i_simd]=0.0;
344  fld.dAs.phi[i_simd]=0.0;
345 #endif
346  }
347 #endif
348 #endif
349  }
350 }
351 
352 
353 template<>
354 template<class Device>
355 KOKKOS_INLINE_FUNCTION void FieldGather<GyroKin>::get_Ah(const ElectricField<Device> fields, const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<GyroKin>& rho, Simd<double>& Ah) const
356 {
357  // Initialize local fields to 0
358 #ifdef EXPLICIT_EM
359  Ah = 0.0;
360 
361  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
362  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
363  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
364  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
365 
366  // Get rho and phi weights
367  FieldWeights<GyroKin> wts(grid, phi, i_simd, rho);
368 #ifdef NEWGYROMATRIX
369  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
370  int irho = wts.rho[0].i;
371 #else
372  int irho = wts.rho.i;
373 #endif
374 
375  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
376  for (int ip = 0; ip<3; ip++){
377  int node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
378  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
379 
380  gather(Ah,i_simd, wp, wts, fields.Ah_rho_ff(node,irho),fields.Ah_rho_ff(node,irho+1));
381  }
382  }
383 #endif
384 }
385 
386 template<>
387 template<class Device>
388 KOKKOS_INLINE_FUNCTION void FieldGather<DriftKin>::get_Ah(const ElectricField<Device> fields, const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<DriftKin>& rho, Simd<double>& Ah) const
389 {
390  // Initialize local fields to 0
391 #ifdef EXPLICIT_EM
392  Ah = 0.0;
393 
394  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
395  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
396  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
397  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
398 
399  // Get rho and phi weights
400  FieldWeights<DriftKin> wts(grid, phi, i_simd, rho);
401 
402  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
403  for (int ip = 0; ip<3; ip++){
404  int node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
405  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
406 
407  gather(Ah,i_simd, wp, wts.phi.w, fields.Ah_phi_ff(wts.phi.i,node));
408  }
409  }
410 #endif
411 }
412 
413 template<>
414 template<class Device>
415 KOKKOS_INLINE_FUNCTION void FieldGather<DriftKin>::get_Ah_cv(const ElectricField<Device> fields, const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, SimdGyroRadius<DriftKin>& rho, Simd<double>& Ah_cv) const
416 {
417  // Initialize local fields to 0
418 #ifdef EXPLICIT_EM
419  Ah_cv = 0.0;
420 
421  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
422  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
423  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
424  double phi= itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
425 
426  // Get rho and phi weights
427  FieldWeights<DriftKin> wts(grid, phi, i_simd, rho);
428 
429  // Loop over the 3 vertices of the grid triangle and add up the contributions from each
430  for (int ip = 0; ip<3; ip++){
431  int node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
432  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
433 
434  gather(Ah_cv,i_simd, wp, wts.phi.w, fields.Ah_cv_phi_ff(wts.phi.i,node));
435  }
436  }
437 #endif
438 }
439 
440 template<KinType PT>
441 template<class Device>
442 KOKKOS_INLINE_FUNCTION void FieldGather<PT>::get_dpot(const ElectricField<Device>& fields, const Grid<Device> &grid, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p, Simd<double>& dpot_out) const {
443 
444  dpot_out.zero();
445 
446  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
447  // If particle is inactive, set itr_work=1, phi=0D0, iphi=0
448  int itr_work = itr[i_simd]>0 ? itr[i_simd] : 1;
449 
450  double phi = itr[i_simd]>0 ? fld_phi[i_simd] : 0.0;
452  FieldWeights<DriftKin> wts(grid, phi, i_simd, rho);
453 
454  for(int ip=0; ip<3; ip++){
455  int node = itr[i_simd]>0 ? (grid.nodes(itr_work - 1).vertex[ip] - 1) : 0; // - 1 since its 0-indexed now
456  double wp = itr[i_simd]>0 ? p[ip][i_simd] : 0.0; // all these trinary ops are probably unnecessary
457 
458 #ifdef XGC1
459 //#ifndef EXPLICIT_EM
460  // for electron -- rho=0 case, optimized.
461  gather(dpot_out,i_simd, wp, wts.phi.w, fields.dpot_ff(node));
462 //#else
463 // dpot_out[i_simd]+=wp*fields.dpot_n0(node);
464 //#endif
465 #elif defined(XGCA)
466  // Simple gather
467  dpot_out[i_simd]+=wp*fields.dpot(node);
468 #endif
469  }
470  }
471 }
472 
473 #endif
Simd< double > r
Definition: globals.hpp:60
Definition: electric_field.hpp:9
Definition: globals.hpp:59
Definition: gyro_radius.hpp:23
KOKKOS_INLINE_FUNCTION void get_Ah_cv(const ElectricField< Device > fields, const Grid< Device > &grid, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, SimdGyroRadius< PT > &rho, Simd< double > &Ah_cv) const
Definition: sml.hpp:8
double V[3][2]
Definition: electric_field.hpp:10
Definition: field_weights.hpp:14
KOKKOS_INLINE_FUNCTION void field_correction(int i_simd, double basis, double gamma_psi, const SimdVector2D &gradpsi, double *rvec, double *zvec)
Definition: field_gather.hpp:13
Definition: electric_field.hpp:35
KOKKOS_INLINE_FUNCTION void get_Ah(const ElectricField< Device > fields, const Grid< Device > &grid, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, SimdGyroRadius< PT > &rho, Simd< double > &Ah) const
Definition: electric_field.hpp:14
Definition: grid.hpp:10
KOKKOS_INLINE_FUNCTION void gather_all_fields(const ElectricField< Device > &fields, const Simulation< Device > &sml, int i_simd, int node, double wp, const double(&rvec)[2], const double(&zvec)[2], T &wts, LocalFields &fld) const
Definition: field_gather.hpp:61
Definition: electric_field.hpp:24
double S[2]
Definition: electric_field.hpp:25
Kokkos::View< FieldXGCa **, Kokkos::LayoutRight, Device > E_rho
Definition: electric_field.hpp:70
LinearWeights rho
Definition: field_weights.hpp:67
Definition: local_fields.hpp:7
r coordinate
Definition: globals.hpp:19
Kokkos::View< int *, Kokkos::LayoutRight, Device > basis
A basis for the guesses?
Definition: grid.hpp:57
Simd< double > z
Definition: globals.hpp:56
double w[2]
Definition: linear_weights.hpp:9
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 VectorFieldPlanes &field) const
double V[2][3]
Definition: electric_field.hpp:15
Definition: field_weights.hpp:65
Simd< double > phi
Definition: globals.hpp:62
Definition: electric_field.hpp:29
Simd< double > z
Definition: globals.hpp:61
int i
Definition: linear_weights.hpp:8
double E[2]
Definition: electric_field.hpp:30
KOKKOS_INLINE_FUNCTION void zero()
Definition: simd.hpp:76
KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector &vec, int i_simd, const SimdVector &B, double Bmag) const
Definition: field_gather.hpp:203
Definition: grid_structs.hpp:7
Simd< double > r
Definition: globals.hpp:55
Definition: electric_field.hpp:19
bool turb_efield
Whether the electric field is turbulent?
Definition: electric_field.hpp:43
KOKKOS_INLINE_FUNCTION void fields_at_point(const ElectricField< Device > fields, const Simulation< Device > &sml, 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: field_gather.hpp:285
Definition: globals.hpp:54
Definition: field_weights.hpp:77
KOKKOS_INLINE_FUNCTION void zero()
Definition: globals.hpp:65
phi coordinate
Definition: globals.hpp:21
KOKKOS_INLINE_FUNCTION void get_dpot(const ElectricField< Device > &fields, const Grid< Device > &grid, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, Simd< double > &dpot_out) const
Definition: field_gather.hpp:442
Definition: gyro_radius.hpp:95
Definition: gyro_radius.hpp:73
Definition: field_gather.hpp:23
z coordinate
Definition: globals.hpp:20
SimdVector E
Definition: local_fields.hpp:8
double V[2][2]
Definition: electric_field.hpp:20
Kokkos::View< double *, Kokkos::LayoutRight, Device > dpot
Definition: electric_field.hpp:71
Kokkos::View< Vertex *, Kokkos::LayoutRight, Device > nodes
Definition: grid.hpp:58