XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
field_aligned_local_fields.hpp
Go to the documentation of this file.
1 #ifndef FIELD_ALIGNED_LOCAL_FIELDS_HPP
2 #define FIELD_ALIGNED_LOCAL_FIELDS_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 #include "field_gather.hpp"
11 
12 
13 /* This class manages how to obtain the local grid field data during the push. It either:
14  * - a) performs a gather operation when the method fields_at_point is used, or
15  * - b) performs a gather operation when initialized, stores those values and
16  * returns them when the fields_at_point method is used.
17  * Option b) is used in XGC1 for ions, because the ion location in an intermediate RK step
18  * might not have a valid associated field since it could be outside of the toroidal sector
19  * for which field information is present.
20  * The generic declaration uses Option a) */
21 template<KinType KT> class FieldAlignedLocalFields{
22  // No members
23 
24  public:
25 
26  // Nothing to do in constructor
27  template<class Device>
28  KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid<Device> &grid, const Simulation<Device> &sml, const Species<Device> &species,
29  const MagneticField<Device> &magnetic_field, const ElectricField<Device> &electric_field, SimdParticles &part, Simd<int>& itr){}
30 
31  /* Gather fields if calculating for electrons */
32  template<class Device>
33  KOKKOS_INLINE_FUNCTION void fields_at_point(const ElectricField<Device> fields, const Simulation<Device> &sml, const Grid<Device> &grid,
34  const SimdVector &B, const SimdVector2D &gradpsi, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p,
35  SimdGyroRadius<KT>& rho, LocalFields& fld) const{
36  FieldGather<KT> fg;
37  fg.fields_at_point(fields, sml, grid, B,gradpsi,fld_phi,itr,p,rho, fld);
38  }
39 };
40 
41 #ifdef XGC1
42 /* If using XGC1, then the class specialization on KinType GyroKin stores the field aligned
43  * local field calculated in the constructor, and converts it to cylindrical coordinates in
44  * fields_at_point. */
45 template<>
47  // Simd vectors are hard-coded to (r, z, phi), but here they represent field-aligned coordinates instead
48  // "phi" --> parallel to B-field
49  SimdVector E;
50 # ifdef EXPLICIT_EM
51  SimdVector dAh;
52  SimdVector dAs;
53 # endif
54 
55  // Members that don't get adjusted but need to be stored
56 #ifdef DELTAF_CONV
57  SimdVector E00;
58  Simd<double> ddpotdt;
59 #endif
60 #ifdef EXPLICIT_EM
61  Simd<double> Ah;
62  Simd<double> As;
63  Simd<double> Epar_em;
64 #endif
65 
66  /* Convert stored electric/magnetic fields to field-aligned coordinates from cylindrical coordinates */
67  KOKKOS_INLINE_FUNCTION void convert_cyl_to_mag(const SimdVector &B, const SimdVector2D &gradpsi, const LocalFields& fld){
68  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
69  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]);
70  double Bp = sqrt(B.r[i_simd]*B.r[i_simd] + B.z[i_simd]*B.z[i_simd]);
71  double grad2 = gradpsi.r[i_simd]*gradpsi.r[i_simd] + gradpsi.z[i_simd]*gradpsi.z[i_simd];
72 
73  // get E-field in magnetic field-aligned basis
74  E.r[i_simd] = (fld.E.r[i_simd] *gradpsi.r[i_simd] + fld.E.z[i_simd] *gradpsi.z[i_simd])/grad2;
75  E.z[i_simd] = (fld.E.r[i_simd] *B.r[i_simd] + fld.E.z[i_simd] *B.z[i_simd] )/Bp;
76  E.phi[i_simd] = (E.z[i_simd]*Bp + fld.E.phi[i_simd]*B.phi[i_simd] )/Bmag; // parallel field
77 
78 # ifdef EXPLICIT_EM
79  // grad(A_h)
80  dAh.r[i_simd] = (fld.dAh.r[i_simd]*gradpsi.r[i_simd] + fld.dAh.z[i_simd]*gradpsi.z[i_simd])/grad2;
81  dAh.z[i_simd] = (fld.dAh.r[i_simd]*B.r[i_simd] + fld.dAh.z[i_simd]*B.z[i_simd] )/Bp;
82  dAh.phi[i_simd] = (dAh.z[i_simd] *Bp + fld.dAh.phi[i_simd]*B.phi[i_simd] )/Bmag;
83 
84  // These variables are for the mixed-variables formulation.
85  // Make sure that they are initialized to 0 if Hamiltonian
86  // formulation is used.
87  dAs.r[i_simd] = (fld.dAs.r[i_simd]*gradpsi.r[i_simd] + fld.dAs.z[i_simd]*gradpsi.z[i_simd])/grad2;
88  dAs.z[i_simd] = (fld.dAs.r[i_simd]*B.r[i_simd] + fld.dAs.z[i_simd]*B.z[i_simd] )/Bp;
89  dAs.phi[i_simd] = (dAs.z[i_simd] *Bp + fld.dAs.phi[i_simd]*B.phi[i_simd] )/Bmag;
90 # endif
91  }
92  }
93 
94  /* Convert stored electric/magnetic fields from field-aligned coordinates to cylindrical coordinates */
95  KOKKOS_INLINE_FUNCTION void convert_mag_to_cyl(const SimdVector &B, const SimdVector2D &gradpsi, LocalFields& fld) const{
96  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
97  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]);
98  double Bp = sqrt(B.r[i_simd]*B.r[i_simd] + B.z[i_simd]*B.z[i_simd]);
99  double Br_norm = B.r[i_simd] / Bp;
100  double Bz_norm = B.z[i_simd] / Bp;
101 
102  fld.E.r[i_simd] = E.r[i_simd]*gradpsi.r[i_simd]+E.z[i_simd]*Br_norm;
103  fld.E.z[i_simd] = E.r[i_simd]*gradpsi.z[i_simd]+E.z[i_simd]*Bz_norm;
104  fld.E.phi[i_simd] = (E.phi[i_simd]*Bmag - fld.E.r[i_simd]*B.r[i_simd] - fld.E.z[i_simd]*B.z[i_simd])/B.phi[i_simd];
105 
106 # ifdef EXPLICIT_EM
107  fld.dAh.r[i_simd] = dAh.r[i_simd]*gradpsi.r[i_simd]+dAh.z[i_simd]*Br_norm;
108  fld.dAh.z[i_simd] = dAh.r[i_simd]*gradpsi.z[i_simd]+dAh.z[i_simd]*Bz_norm;
109  fld.dAh.phi[i_simd] = (dAh.phi[i_simd]*Bmag - fld.dAh.r[i_simd]*B.r[i_simd] - fld.dAh.z[i_simd]*B.z[i_simd])/B.phi[i_simd];
110 
111  // The following are for the mixed-variable formulation
112  // Make sure that no all variables are initialized
113  fld.dAs.r[i_simd] = dAs.r[i_simd]*gradpsi.r[i_simd]+dAs.z[i_simd]*Br_norm;
114  fld.dAs.z[i_simd] = dAs.r[i_simd]*gradpsi.z[i_simd]+dAs.z[i_simd]*Bz_norm;
115  fld.dAs.phi[i_simd] = (dAs.phi[i_simd]*Bmag - fld.dAs.r[i_simd]*B.r[i_simd] - fld.dAs.z[i_simd]*B.z[i_simd])/B.phi[i_simd];
116 # endif
117  }
118  }
119 
120  public:
121 
122  /* Gather electric/magnetic fields and convert to field-aligned coordinates for storage */
123  template<class Device>
124  KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid<Device> &grid, const Simulation<Device> &sml, const Species<Device> &species,
125  const MagneticField<Device> &magnetic_field, const ElectricField<Device> &electric_field, SimdParticles &part, Simd<int>& itr){
126 
127  // Would be cleaner to just pass r and z instead of consolidating
128  SimdVector2D x;
129  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
130  x.r[i_simd] = part.ph.r[i_simd];
131  x.z[i_simd] = part.ph.z[i_simd];
132  }
133 
134  // Declare arrays returned from B-field
135  SimdVector B;
136  SimdVector jacb[3];
137  Simd<double> psi;
138  SimdVector2D gradpsi;
139  SimdVector tdb;
140  Simd<bool> rz_outside;
141 
142  // obtain B-field information
143  magnetic_field.field(part.ph.r, part.ph.z, B, jacb, psi, gradpsi, tdb, rz_outside);
144  remove_particles(part.gid,rz_outside);
145 
146  // Get toroidal angle
147  Simd<double> fld_phi;
148  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
149  fld_phi[i_simd] = part.ph.phi[i_simd];
150  }
151  grid.wedge_modulo_phi(fld_phi);
152 
153  // find triangle
154  SimdGridVec p;
155  SimdVector2D xff;
156  grid.charge_search_index(magnetic_field, x, part.ph.phi,psi,xff,itr,p);
157 
158  // calculate rho
159  SimdGyroRadius<GyroKin> rho(grid, magnetic_field, species, x, fld_phi, part.ct.mu, itr);
160 
161  // get electric field
162  LocalFields fld;
164  fg.fields_at_point(electric_field, sml, grid, B,gradpsi,fld_phi,itr,p,rho, fld);
165 
166  // Convert to field-aligned coordinates and store
167  convert_cyl_to_mag(B, gradpsi, fld);
168 
169  // These get stored without any conversion; Consolidate field_aligned_local_fields with local_fields objects to avoid this copy
170 #ifdef DELTAF_CONV
171  E00.r=fld.E00.r; // Should E00 also be field-aligned?
172  E00.z=fld.E00.z;
173  E00.phi=fld.E00.phi;
174  ddpotdt = fld.ddpotdt;
175 #endif
176 #ifdef EXPLICIT_EM
177  Ah = fld.Ah;
178  As = fld.As;
179  Epar_em = fld.Epar_em;
180 #endif
181  }
182 
183  /* Return the stored field values for XGC1 ions */
184  template<class Device>
185  KOKKOS_INLINE_FUNCTION void fields_at_point(const ElectricField<Device> fields, const Simulation<Device> &sml, const Grid<Device> &grid,
186  const SimdVector &B, const SimdVector2D &gradpsi, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p,
187  SimdGyroRadius<GyroKin>& rho, LocalFields& fld) const{
188  convert_mag_to_cyl(B, gradpsi, fld);
189 
190  // These get stored without any conversion; Consolidate field_aligned_local_fields with local_fields objects to avoid this copy
191 #ifdef DELTAF_CONV
192  fld.E00.r=E00.r; // Should E00 also be field-aligned?
193  fld.E00.z=E00.z;
194  fld.E00.phi=E00.phi;
195  fld.ddpotdt = ddpotdt;
196 #endif
197 #ifdef EXPLICIT_EM
198  fld.Ah = Ah;
199  fld.As = As;
200  fld.Epar_em = Epar_em;
201 #endif
202  }
203 };
204 #endif
205 
206 
207 #endif
Simd< double > r
Definition: globals.hpp:60
Definition: globals.hpp:59
Definition: gyro_radius.hpp:23
Definition: sml.hpp:8
KOKKOS_INLINE_FUNCTION void field(const Simd< double > &fld_r, const Simd< double > &fld_z, SimdVector &bvec, SimdVector(&jacb)[3], Simd< double > &psivec, SimdVector2D &gradpsi, SimdVector &tdb, Simd< bool > &rz_outside) const
Definition: magnetic_field.tpp:78
Definition: magnetic_field.hpp:9
Definition: electric_field.hpp:35
Definition: grid.hpp:10
Definition: local_fields.hpp:7
Simd< double > r
Definition: particles.hpp:15
Simd< double > z
Definition: globals.hpp:56
Definition: globals.hpp:15
SimdPhase ph
Definition: particles.hpp:32
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< KT > &rho, LocalFields &fld) const
Definition: field_aligned_local_fields.hpp:33
KOKKOS_INLINE_FUNCTION void wedge_modulo_phi(Simd< double > &phi_mod) const
Definition: grid.tpp:554
Simd< double > phi
Definition: globals.hpp:62
Definition: particles.hpp:31
Simd< double > z
Definition: globals.hpp:61
KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid< Device > &grid, const Simulation< Device > &sml, const Species< Device > &species, const MagneticField< Device > &magnetic_field, const ElectricField< Device > &electric_field, SimdParticles &part, Simd< int > &itr)
Definition: field_aligned_local_fields.hpp:28
Definition: grid_structs.hpp:7
Simd< double > z
Definition: particles.hpp:16
Simd< double > r
Definition: globals.hpp:55
Definition: field_aligned_local_fields.hpp:21
Simd< long long int > gid
Definition: particles.hpp:34
Simd< double > phi
Definition: particles.hpp:17
Definition: magnetic_field.F90:1
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: simd.hpp:18
Definition: globals.hpp:54
SimdConstants ct
Definition: particles.hpp:33
Definition: species.hpp:13
Definition: gyro_radius.hpp:73
Definition: field_gather.hpp:23
Simd< double > mu
Definition: particles.hpp:25
SimdVector E
Definition: local_fields.hpp:8
KOKKOS_INLINE_FUNCTION void remove_particles(Simd< long long int > &gid, const Simd< bool > &deactivate)
Definition: particles.tpp:4
KOKKOS_INLINE_FUNCTION void charge_search_index(const MagneticField< Device > &magnetic_field, const SimdVector2D &x, const Simd< double > &phi, const Simd< double > &psi, SimdVector2D &xff, Simd< int > &itr, SimdGridVec &p) const
Definition: grid.tpp:505