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 "push_controls.hpp"
6 #include "grid.hpp"
7 #include "field_weights.hpp"
8 #include "gyro_radius.hpp"
9 #include "electric_field.hpp"
10 #include "grid_field_pack.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, PhiInterpType PIT>
23  // No members
24 
25  public:
26 
27  // Nothing to do in constructor
28  template<class Device>
29  KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid<Device> &grid, const PushControls &push_controls, const Species<Device> &species,
31 
32  /* Gather fields if calculating for electrons */
33  template<class Device>
34  KOKKOS_INLINE_FUNCTION void fields_at_point(const GridFieldPack<Device, PIT> fields, const PushControls &push_controls, const Grid<Device> &grid,
35  const SimdVector &B, const SimdVector2D &gradpsi, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p,
36  SimdGyroRadius<KT>& rho, LocalFields& fld) const{
37  fields.fields_at_point(push_controls, grid, B,gradpsi,fld_phi,itr,p,rho, fld);
38  }
39 };
40 
41 /* If using XGC1, then the class specialization on KinType GyroKin stores the field aligned
42  * local field calculated in the constructor, and converts it to cylindrical coordinates in
43  * fields_at_point. */
44 template<>
46  // Simd vectors are hard-coded to (r, z, phi), but here they represent field-aligned coordinates instead
47  // "phi" --> parallel to B-field
49 # ifdef EXPLICIT_EM
50  SimdVector dAh;
51  SimdVector dAs;
52 # endif
53 
54  // Members that don't get adjusted but need to be stored
55 #ifdef DELTAF_CONV
56  SimdVector E00;
57  Simd<double> ddpotdt;
58 #endif
59 #ifdef EXPLICIT_EM
60  Simd<double> Ah;
61  Simd<double> As;
62  Simd<double> Epar_em;
63 #endif
64 
65  /* Convert stored electric/magnetic fields to field-aligned coordinates from cylindrical coordinates */
66  KOKKOS_INLINE_FUNCTION void convert_cyl_to_mag(const SimdVector &B, const SimdVector2D &gradpsi, const LocalFields& fld){
67  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
68  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]);
69  double Bp = sqrt(B.r[i_simd]*B.r[i_simd] + B.z[i_simd]*B.z[i_simd]);
70  double grad2 = gradpsi.r[i_simd]*gradpsi.r[i_simd] + gradpsi.z[i_simd]*gradpsi.z[i_simd];
71 
72  // get E-field in magnetic field-aligned basis
73  E.r[i_simd] = (fld.E.r[i_simd] *gradpsi.r[i_simd] + fld.E.z[i_simd] *gradpsi.z[i_simd])/grad2;
74  E.z[i_simd] = (fld.E.r[i_simd] *B.r[i_simd] + fld.E.z[i_simd] *B.z[i_simd] )/Bp;
75  E.phi[i_simd] = (E.z[i_simd]*Bp + fld.E.phi[i_simd]*B.phi[i_simd] )/Bmag; // parallel field
76 
77 # ifdef EXPLICIT_EM
78  // grad(A_h)
79  dAh.r[i_simd] = (fld.dAh.r[i_simd]*gradpsi.r[i_simd] + fld.dAh.z[i_simd]*gradpsi.z[i_simd])/grad2;
80  dAh.z[i_simd] = (fld.dAh.r[i_simd]*B.r[i_simd] + fld.dAh.z[i_simd]*B.z[i_simd] )/Bp;
81  dAh.phi[i_simd] = (dAh.z[i_simd] *Bp + fld.dAh.phi[i_simd]*B.phi[i_simd] )/Bmag;
82 
83  // These variables are for the mixed-variables formulation.
84  // Make sure that they are initialized to 0 if Hamiltonian
85  // formulation is used.
86  dAs.r[i_simd] = (fld.dAs.r[i_simd]*gradpsi.r[i_simd] + fld.dAs.z[i_simd]*gradpsi.z[i_simd])/grad2;
87  dAs.z[i_simd] = (fld.dAs.r[i_simd]*B.r[i_simd] + fld.dAs.z[i_simd]*B.z[i_simd] )/Bp;
88  dAs.phi[i_simd] = (dAs.z[i_simd] *Bp + fld.dAs.phi[i_simd]*B.phi[i_simd] )/Bmag;
89 # endif
90  }
91  }
92 
93  /* Convert stored electric/magnetic fields from field-aligned coordinates to cylindrical coordinates */
94  KOKKOS_INLINE_FUNCTION void convert_mag_to_cyl(const SimdVector &B, const SimdVector2D &gradpsi, LocalFields& fld) const{
95  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
96  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]);
97  double Bp = sqrt(B.r[i_simd]*B.r[i_simd] + B.z[i_simd]*B.z[i_simd]);
98  double Br_norm = B.r[i_simd] / Bp;
99  double Bz_norm = B.z[i_simd] / Bp;
100 
101  fld.E.r[i_simd] = E.r[i_simd]*gradpsi.r[i_simd]+E.z[i_simd]*Br_norm;
102  fld.E.z[i_simd] = E.r[i_simd]*gradpsi.z[i_simd]+E.z[i_simd]*Bz_norm;
103  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];
104 
105 # ifdef EXPLICIT_EM
106  fld.dAh.r[i_simd] = dAh.r[i_simd]*gradpsi.r[i_simd]+dAh.z[i_simd]*Br_norm;
107  fld.dAh.z[i_simd] = dAh.r[i_simd]*gradpsi.z[i_simd]+dAh.z[i_simd]*Bz_norm;
108  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];
109 
110  // The following are for the mixed-variable formulation
111  // Make sure that no all variables are initialized
112  fld.dAs.r[i_simd] = dAs.r[i_simd]*gradpsi.r[i_simd]+dAs.z[i_simd]*Br_norm;
113  fld.dAs.z[i_simd] = dAs.r[i_simd]*gradpsi.z[i_simd]+dAs.z[i_simd]*Bz_norm;
114  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];
115 # endif
116  }
117  }
118 
119  public:
120 
121  /* Gather electric/magnetic fields and convert to field-aligned coordinates for storage */
122  template<class Device>
123  KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid<Device> &grid, const PushControls &push_controls, const Species<Device> &species,
125  // Declare arrays returned from B-field
126  SimdVector B;
127  SimdVector jacb[3];
128  Simd<double> psi;
129  SimdVector2D gradpsi;
130  SimdVector tdb;
131  Simd<bool> rz_outside;
132 
133  // obtain B-field information
134  magnetic_field.field(part.ph.x(), B, jacb, psi, gradpsi, tdb, rz_outside);
135  remove_particles(part.gid,rz_outside);
136 
137  // Get toroidal angle
138  Simd<double> fld_phi = part.ph.phi;
139  grid.wedge_modulo_phi(fld_phi);
140 
141  // Get the current field following grid triangle and p weights
142  SimdGridVec p;
143  grid.charge_search_index(magnetic_field, part.ph.v(), psi, itr, p);
144 
145  // calculate rho
146  SimdGyroRadius<GyroKin> rho(grid, magnetic_field, species, part.ph.x(), fld_phi, part.ct.mu, itr);
147 
148  // get electric field
149  LocalFields fld;
150  gfpack.fields_at_point(push_controls, grid, B,gradpsi,fld_phi,itr,p,rho, fld);
151 
152  // Convert to field-aligned coordinates and store
153  convert_cyl_to_mag(B, gradpsi, fld);
154 
155  // These get stored without any conversion; Consolidate field_aligned_local_fields with local_fields objects to avoid this copy
156 #ifdef DELTAF_CONV
157  E00.r=fld.E00.r; // Should E00 also be field-aligned?
158  E00.z=fld.E00.z;
159  E00.phi=fld.E00.phi;
160  ddpotdt = fld.ddpotdt;
161 #endif
162 #ifdef EXPLICIT_EM
163  Ah = fld.Ah;
164  As = fld.As;
165  Epar_em = fld.Epar_em;
166 #endif
167  }
168 
169  /* Return the stored field values for XGC1 ions */
170  template<class Device>
171  KOKKOS_INLINE_FUNCTION void fields_at_point(const GridFieldPack<Device, PhiInterpType::Planes> fields, const PushControls &push_controls, const Grid<Device> &grid,
172  const SimdVector &B, const SimdVector2D &gradpsi, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p,
173  SimdGyroRadius<GyroKin>& rho, LocalFields& fld) const{
174  convert_mag_to_cyl(B, gradpsi, fld);
175 
176  // These get stored without any conversion; Consolidate field_aligned_local_fields with local_fields objects to avoid this copy
177 #ifdef DELTAF_CONV
178  fld.E00.r=E00.r; // Should E00 also be field-aligned?
179  fld.E00.z=E00.z;
180  fld.E00.phi=E00.phi;
181  fld.ddpotdt = ddpotdt;
182 #endif
183 #ifdef EXPLICIT_EM
184  fld.Ah = Ah;
185  fld.As = As;
186  fld.Epar_em = Epar_em;
187 #endif
188  }
189 };
190 
191 
192 #endif
Simd< double > r
Definition: simd.hpp:145
Definition: simd.hpp:144
Definition: gyro_radius.hpp:24
Definition: push_controls.hpp:8
Definition: magnetic_field.hpp:9
Definition: grid.hpp:10
Definition: grid_field_pack.hpp:239
KOKKOS_INLINE_FUNCTION void fields_at_point(const GridFieldPack< Device, PIT > fields, const PushControls &push_controls, 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:34
Definition: local_fields.hpp:7
PhiInterpType
Definition: globals.hpp:88
KOKKOS_INLINE_FUNCTION SimdVector2D & x()
Definition: particles.hpp:27
KOKKOS_INLINE_FUNCTION SimdVector & v()
Definition: particles.hpp:39
KOKKOS_INLINE_FUNCTION void convert_cyl_to_mag(const SimdVector &B, const SimdVector2D &gradpsi, const LocalFields &fld)
Definition: field_aligned_local_fields.hpp:66
Simd< double > z
Definition: simd.hpp:141
Definition: globals.hpp:83
SimdPhase ph
Definition: particles.hpp:59
KOKKOS_INLINE_FUNCTION void wedge_modulo_phi(Simd< double > &phi_mod) const
Definition: grid.tpp:611
Simd< double > phi
Definition: simd.hpp:147
KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid< Device > &grid, const PushControls &push_controls, const Species< Device > &species, const MagneticField< Device > &magnetic_field, const GridFieldPack< Device, PhiInterpType::Planes > &gfpack, SimdParticles &part, Simd< int > &itr)
Definition: field_aligned_local_fields.hpp:123
Definition: particles.hpp:58
Simd< double > z
Definition: simd.hpp:146
Definition: grid_structs.hpp:7
Simd< double > r
Definition: simd.hpp:140
Definition: field_aligned_local_fields.hpp:22
Simd< long long int > gid
Definition: particles.hpp:61
Simd< double > phi
Definition: particles.hpp:20
KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid< Device > &grid, const PushControls &push_controls, const Species< Device > &species, const MagneticField< Device > &magnetic_field, const GridFieldPack< Device, PIT > &gfpack, SimdParticles &part, Simd< int > &itr)
Definition: field_aligned_local_fields.hpp:29
Definition: magnetic_field.F90:1
SimdVector E
Definition: field_aligned_local_fields.hpp:48
Definition: simd.hpp:18
Definition: simd.hpp:139
KOKKOS_INLINE_FUNCTION void field(const SimdVector2D &x, SimdVector &bvec, SimdVector(&jacb)[3], Simd< double > &psivec, SimdVector2D &gradpsi, SimdVector &tdb, Simd< bool > &rz_outside) const
Definition: magnetic_field.tpp:78
SimdConstants ct
Definition: particles.hpp:60
KOKKOS_INLINE_FUNCTION void convert_mag_to_cyl(const SimdVector &B, const SimdVector2D &gradpsi, LocalFields &fld) const
Definition: field_aligned_local_fields.hpp:94
Definition: species.hpp:66
Definition: gyro_radius.hpp:74
Simd< double > mu
Definition: particles.hpp:52
KOKKOS_INLINE_FUNCTION void fields_at_point(const GridFieldPack< Device, PhiInterpType::Planes > fields, const PushControls &push_controls, const Grid< Device > &grid, const SimdVector &B, const SimdVector2D &gradpsi, const Simd< double > &fld_phi, const Simd< int > &itr, const SimdGridVec &p, SimdGyroRadius< GyroKin > &rho, LocalFields &fld) const
Definition: field_aligned_local_fields.hpp:171
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:529
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 Simd< int > &itr, const SimdGridVec &p, SimdGyroRadius< PT > &rho, LocalFields &fld) const
Definition: grid_field_pack.hpp:398