XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
get_drift_velocity.hpp
Go to the documentation of this file.
1 #ifndef GET_DRIFT_VELOCITY_HPP
2 #define GET_DRIFT_VELOCITY_HPP
3 
4 #include "push_physics.hpp"
5 
6 // Modified from fields_at_point in grid_field_pack. Should consolidate and use GridWeights<Order::Zero>
7 template<typename GFPT>
8 KOKKOS_INLINE_FUNCTION void fields_at_node(const GridFieldPack<DeviceType, GFPT>& gfpack, const PushControls &push_controls, const Grid<DeviceType> &grid, const SimdVector &B, const SimdVector2D &gradpsi, const Simd<int>& global_inds, const SimdGridWeights<Order::One, PIT_GLOBAL> grid_wts, const SimdGyroWeights<KinType::GyroKin>& rho_wts, LocalFields<GFPT>& fld){
9  // Initialize local fields to 0
10  fld.template get<Label::E>().zero();
11 #ifdef EXPLICIT_EM
12  fld.template get<Label::dAh>().zero();
13  fld.template get<Label::dAs>().zero();
14  fld.template get<Label::Ah>().zero();
15  fld.template get<Label::As>().zero();
16  fld.template get<Label::Epar_em>().zero();
17 #endif
18 
19  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
20  // Field basis correction
21  FieldCorrection field_correction(gradpsi, i_simd);
22 
23  constexpr double wp = 1.0; // No triangle weights since this is on a vertex
24 
25  field_correction.set(i_simd, grid.uses_rz_basis(global_inds[i_simd]), gradpsi);
26 
27  gfpack.gather_all_fields(push_controls, i_simd, global_inds[i_simd], wp, field_correction, grid_wts, rho_wts, fld);
28 
29 #ifdef XGC1
30  // Correct phi component of E-field
31  double Bmag = B.magnitude(i_simd);
32 
33  // Get phi from para
34  gfpack.phi_from_para(fld.template get<Label::E>(), i_simd, B, Bmag);
35 #ifdef EXPLICIT_EM
36  gfpack.phi_from_para(fld.template get<Label::dAh>(), i_simd, B, Bmag);
37  gfpack.phi_from_para(fld.template get<Label::dAs>(), i_simd, B, Bmag);
38 #endif
39 #endif
40  }
41 }
42 
43 // Modified from fields_at_point in grid_field_pack. Should consolidate and use GridWeights<Order::Zero>
45  // Initialize local fields to 0
46  pot_out = 0.0;
47 #ifdef XGC1
48  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
49  int irho = rho_wts.irho(i_simd);
50  constexpr double wp = 1.0; // No triangle weights since this is on a vertex
51  pot(global_inds[i_simd], irho).gather(pot_out,i_simd, wp, rho_wts, grid_wts.phi_wts, pot(global_inds[i_simd], irho+1));
52  }
53 #endif
54 }
55 
56 // Parallel drift
57 // v_parallel = D*(B_star*dH/dv_p *1/(m*B) + v_p*dB_RMP)
58 // But we move the m/c*u_para*curl(b) part to the magnetic drift term.
59 // The remaining terms form the actual parallel motion along
60 // the perturbed fieldn --> slightly different B_star
61 template<typename GFPT>
62 KOKKOS_INLINE_FUNCTION void get_v_par_drift(double D, const SimdVector& bfield, const SimdVector& tdb, const double (&curl_nb)[3], const LocalFields<GFPT>& fld, double over_B, double cmrho, double gradpsi_r_norm, double gradpsi_z_norm, double (&v_par)[3], int i_simd){
63  // B_star = B+A_s*curl(b) + grad(A)xB
64  double b_star[3];
65  b_star[PIR] = bfield.r[i_simd];
66  b_star[PIZ] = bfield.z[i_simd];
67  b_star[PIP] = bfield.phi[i_simd];
68 #ifdef EXPLICIT_EM
69  const auto& As = fld.template get<Label::As>();
70  const auto& dAs = fld.template get<Label::dAs>();
71  constexpr int NEG_DA = -1; // dAs needs sign reversal!
72  b_star[PIR] += As[i_simd]*curl_nb[PIR] + NEG_DA*(dAs.phi[i_simd]*bfield.z[i_simd] - dAs.z[i_simd]*bfield.phi[i_simd])*over_B;
73  b_star[PIZ] += As[i_simd]*curl_nb[PIZ] + NEG_DA*(dAs.r[i_simd]*bfield.phi[i_simd] - dAs.phi[i_simd]*bfield.r[i_simd])*over_B;
74  b_star[PIP] += As[i_simd]*curl_nb[PIP] + NEG_DA*(dAs.z[i_simd]*bfield.r[i_simd] - dAs.r[i_simd]*bfield.z[i_simd])*over_B;
75 #else
76  b_star[PIR] += tdb.r[i_simd];
77  b_star[PIZ] += tdb.z[i_simd];
78 #endif
79 
80  // In the EM formalism, the parallel flow is v_para*B_star, i.e., he parallel flow
81  // contribution to the radial flux can come from the vr and dvr --->
82  // Subtract the curl(B) contribution, which is accounted for in v_mag
83  double prefac = D*cmrho;
84  v_par[PIR] = prefac*b_star[PIR];
85  v_par[PIZ] = prefac*b_star[PIZ];
86  v_par[PIP] = prefac*b_star[PIP];
87 }
88 
89 // The magnetic cross-field drifts
90 // Remove the (numerical flux-surface average of the magnetic drifts.
91 // This approximation accepts a small error because the velocity dependence
92 // in "D" is not taken into accound in v_gradb_avg and v_curv_avg.
93 //
94 // v_drift = D*(FxB) = D*(Bxgrad(H))/B^2 !Use all terms except E in grad(H)
95 // f = f_mr + f_pert = f_mr + f_es + f_em
96 // f_mr --> Mirror force
97 // f_es --> Electric field
98 // f_em --> Electromagnetic forces
99 // relevant force f = f_mr + f_em
100 // = -mu/c * grad(B) + (v_para-c/m*Ah)*grad(Ah)
101 // -mu/c*[grad(B)xB]/B^2 --> Apply cancellation law
102 template<typename GFPT>
103 KOKKOS_INLINE_FUNCTION void get_v_mag_drift(double D, double cmrho, double mu, double u, double rho, double c_m, double charge, double inv_r, const double (&curl_B)[3], const double (&grad_B)[3], const LocalFields<GFPT>& fld, const SimdVector& bfield, double B, double over_B2, double gradpsi_r_norm, double gradpsi_z_norm, const Simd<double>& v_curv_rad_fsa, const Simd<double>& v_grad_B_rad_fsa, double (&v_mag)[3], int i_simd){
104  // Convert to r,z
105  double v_curv_r_fsa, v_curv_z_fsa, v_grad_B_r_fsa, v_grad_B_z_fsa;
106  psitheta_2_rz(gradpsi_r_norm, gradpsi_z_norm, v_curv_rad_fsa[i_simd], 0.0, v_curv_r_fsa, v_curv_z_fsa);
107  psitheta_2_rz(gradpsi_r_norm, gradpsi_z_norm, v_grad_B_rad_fsa[i_simd], 0.0, v_grad_B_r_fsa, v_grad_B_z_fsa);
108 
109  double over_c = 1.0/charge;
110 
111  // Grad B drift is Bxgrad(|B|)
112  double grad_B_drift[3];
113  grad_B_drift[PIR] = bfield.phi[i_simd] * grad_B[PIZ] - bfield.z[i_simd] * grad_B[PIP]*inv_r;
114  grad_B_drift[PIZ] = bfield.r[i_simd] * grad_B[PIP]*inv_r - bfield.phi[i_simd] * grad_B[PIR];
115  grad_B_drift[PIP] = bfield.z[i_simd] * grad_B[PIR] - bfield.r[i_simd] * grad_B[PIZ];
116 
117 #ifdef EXPLICIT_EM
118  double over_B = 1.0/B;
119  const auto& Ah = fld.template get<Label::Ah>();
120  const auto& dAh = fld.template get<Label::dAh>();
121  constexpr int NEG_DA = -1; // dAh needs sign reversal!
122  double fr = -mu*over_c*grad_B[PIR] + NEG_DA*u*dAh.r[i_simd];
123  double fz = -mu*over_c*grad_B[PIZ] + NEG_DA*u*dAh.z[i_simd];
124  double fp = -mu*over_c*grad_B[PIP] + NEG_DA*u*dAh.phi[i_simd];
125 #ifdef EXPLICIT_EM_NONLIN_TERMS
126  fr -= NEG_DA*c_m*Ah[i_simd]*dAh.r[i_simd];
127  fz -= NEG_DA*c_m*Ah[i_simd]*dAh.z[i_simd];
128  fp -= NEG_DA*c_m*Ah[i_simd]*dAh.phi[i_simd];
129 #endif
130 
131  // Magnetic drifts in (R,phi,Z) coordinates:
132  v_mag[PIR] = D*(bfield.z[i_simd]*fp - bfield.phi[i_simd]*fz)*over_B2;
133  v_mag[PIZ] = D*(bfield.phi[i_simd]*fr - bfield.r[i_simd]*fp)*over_B2;
134  v_mag[PIP] = D*(bfield.r[i_simd]*fz - bfield.z[i_simd]*fr)*over_B2;
135 
136  // Apply cancellation law of Bxgrad(B)
137  // Also add contribution from curl(b) (removed from B*dH/du)
138  v_mag[PIR] -= D*over_c*over_B2*mu*v_grad_B_r_fsa;
139  v_mag[PIZ] -= D*over_c*over_B2*mu*v_grad_B_z_fsa;
140 
141  v_mag[PIR] += D*over_B*u*rho*(over_B*(grad_B_drift[PIR] - v_grad_B_r_fsa) + (B*curl_B[PIR] - v_curv_r_fsa));
142  v_mag[PIZ] += D*over_B*u*rho*(over_B*(grad_B_drift[PIZ] - v_grad_B_z_fsa) + (B*curl_B[PIZ] - v_curv_z_fsa));
143  v_mag[PIP] += D*over_B*u*rho*(over_B* grad_B_drift[PIP] + B*curl_B[PIP]);
144 #else
145  double cmrho2 = cmrho*rho;
146  double murho2B = (mu + charge*cmrho2*B);
147  double murho2B_c = murho2B*over_c;
148 
149  // The magnetic drifts (curvature + grad_b)
150  // Remove the (numerical flux-surface average of the magnetic drifts.
151  // This approximation accepts a small error because the velocity dependence
152  // in "D" is not taken into accound in v_gradb_avg and v_curv_avg.
153  v_mag[PIR] = D*((grad_B_drift[PIR] - v_grad_B_r_fsa)*murho2B_c*over_B2 + (curl_B[PIR] - v_curv_r_fsa)*cmrho2);
154  v_mag[PIZ] = D*((grad_B_drift[PIZ] - v_grad_B_z_fsa)*murho2B_c*over_B2 + (curl_B[PIZ] - v_curv_z_fsa)*cmrho2);
155  v_mag[PIP] = D*( grad_B_drift[PIP] *murho2B_c*over_B2 + curl_B[PIP] *cmrho2);
156 #endif
157 }
158 
159 
160 //> Evaluates the drift velocities at given inode, imu, ivp
161 // The basic magnetic field quantities necessary for the calculation of the
162 // magnetic drifts are stored in the grid data structure.
163 // To save memory, this could also be shifted to f0_module.
164 //
165 // CAUTION: Although this routine uses the full expression for dR/dt,
166 // it is correct only when using the pullback method and when
167 // in addition, A_h=0, i.e., directly after the pullback operation
168 // because it uses the assumption that u_para=v_para
169 //
170 template<typename GFPT>
171 KOKKOS_INLINE_FUNCTION void get_drift_velocity(const Species<DeviceType>& species, const Simd<double>& vp, const Simd<double>& v_curv_rad_fsa, const Simd<double>& v_grad_B_rad_fsa, const Simd<double>& mu, const SimdVector& bfield, const SimdVector (&jacb)[3], const SimdVector& tdb, const Simd<double>& B_mag, const Simd<double>& inv_r, const Simd<double>& gradpsi_r_norm, const Simd<double>& gradpsi_z_norm, const LocalFields<GFPT>& fld, SimdVector& v_mag, SimdVector& v_exb, SimdVector& v_par){
172 
173  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
174  const double& B = B_mag[i_simd];
175  double rho = vp[i_simd]/(B*species.c_m);
176 
177  double over_B = 1.0/B;
178  double over_B2 = over_B*over_B;
179  double cmrho = species.c_m*rho;
180 
181  // For clarity, distinguish v_para and u_para
182  double u = vp[i_simd];
183 #ifdef EXPLICIT_EM
184  u += species.c_m*fld.template get<Label::Ah>()[i_simd];
185 #endif
186 
187  double grad_B[3];
188  get_grad_B(bfield, jacb, over_B, grad_B, i_simd);
189 
190  double curl_B[3];
191  get_curl_B(bfield, jacb, inv_r[i_simd], curl_B, i_simd);
192 
193  double curl_nb[3]; // [ curl of normalized vec B ] = curl vec B / B + vec B X grad B / B^2
194  get_curl_nb(bfield, grad_B, over_B, over_B2, curl_B, inv_r[i_simd], curl_nb, i_simd);
195 
196  // Velocity space Jacobian factor
197 #ifdef EXPLICIT_EM
198  const double D = get_vspace_jacb_factor(rho + fld.template get<Label::As>()[i_simd] * over_B, over_B2, bfield, jacb, inv_r[i_simd], i_simd);
199 #else
200  const double D = get_vspace_jacb_factor(rho, over_B2, bfield, jacb, inv_r[i_simd], i_simd);
201 #endif
202 
203  // The ExB drift
204  constexpr double drift_on = 1.0;
205  double yp[3];
206  get_ExB_drift(D, drift_on, bfield, fld.template get<Label::E>(), i_simd, over_B2, yp);
207  v_exb.r[i_simd] = yp[PIR];
208  v_exb.z[i_simd] = yp[PIZ];
209  v_exb.phi[i_simd] = yp[PIP];
210 
211  // Magnetic drifts
212  get_v_mag_drift(D, cmrho, mu[i_simd], u, rho, species.c_m, species.charge, inv_r[i_simd], curl_B, grad_B, fld, bfield, B, over_B2, gradpsi_r_norm[i_simd], gradpsi_z_norm[i_simd], v_curv_rad_fsa, v_grad_B_rad_fsa, yp, i_simd);
213  v_mag.r[i_simd] = yp[PIR];
214  v_mag.z[i_simd] = yp[PIZ];
215  v_mag.phi[i_simd] = yp[PIP];
216 
217  // Parallel drift
218  get_v_par_drift(D, bfield, tdb, curl_nb, fld, over_B, cmrho, gradpsi_r_norm[i_simd], gradpsi_z_norm[i_simd], yp, i_simd);
219  v_par.r[i_simd] = yp[PIR];
220  v_par.z[i_simd] = yp[PIZ];
221  v_par.phi[i_simd] = yp[PIP];
222  }
223 }
224 
225 #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_grad_B(const SimdVector &bfield, const SimdVector(&jacb)[3], double over_B, double(&grad_B)[3], int i_simd)
Definition: push_physics.tpp:155
KOKKOS_INLINE_FUNCTION void get_drift_velocity(const Species< DeviceType > &species, const Simd< double > &vp, const Simd< double > &v_curv_rad_fsa, const Simd< double > &v_grad_B_rad_fsa, const Simd< double > &mu, const SimdVector &bfield, const SimdVector(&jacb)[3], const SimdVector &tdb, const Simd< double > &B_mag, const Simd< double > &inv_r, const Simd< double > &gradpsi_r_norm, const Simd< double > &gradpsi_z_norm, const LocalFields< GFPT > &fld, SimdVector &v_mag, SimdVector &v_exb, SimdVector &v_par)
Definition: get_drift_velocity.hpp:171
double c_m
c/m
Definition: species.hpp:86
KOKKOS_INLINE_FUNCTION void get_v_mag_drift(double D, double cmrho, double mu, double u, double rho, double c_m, double charge, double inv_r, const double(&curl_B)[3], const double(&grad_B)[3], const LocalFields< GFPT > &fld, const SimdVector &bfield, double B, double over_B2, double gradpsi_r_norm, double gradpsi_z_norm, const Simd< double > &v_curv_rad_fsa, const Simd< double > &v_grad_B_rad_fsa, double(&v_mag)[3], int i_simd)
Definition: get_drift_velocity.hpp:103
Definition: grid_weights.hpp:47
KOKKOS_INLINE_FUNCTION double get_vspace_jacb_factor(double rho_mod, double over_B2, const SimdVector &bfield, const SimdVector(&jacb)[3], double inv_r, int i_simd)
Definition: push_physics.tpp:137
Definition: push_controls.hpp:9
Definition: grid_field_pack.hpp:24
KOKKOS_INLINE_FUNCTION void get_curl_B(const SimdVector &bfield, const SimdVector(&jacb)[3], double inv_r, double(&curl_B)[3], int i_simd)
Definition: push_physics.tpp:102
KOKKOS_INLINE_FUNCTION int uses_rz_basis(const int inode) const
Definition: grid.tpp:252
Definition: grid_field.hpp:22
Definition: gyro_radius.hpp:24
r coordinate
Definition: globals.hpp:201
Definition: field.hpp:23
double charge
Particle charge.
Definition: species.hpp:84
KOKKOS_INLINE_FUNCTION void pot_rho_at_node(const GridField< DeviceType, VarType::Scalar, PIT_GLOBAL, TorType::OnePlane, KinType::GyroKin > &pot, const Grid< DeviceType > &grid, const Simd< int > &global_inds, const SimdGridWeights< Order::One, PIT_GLOBAL > grid_wts, const SimdGyroWeights< KinType::GyroKin > &rho_wts, Simd< double > &pot_out)
Definition: get_drift_velocity.hpp:44
KOKKOS_INLINE_FUNCTION void psitheta_2_rz(double dpsidr, double dpsidz, double x_rad, double x_theta, double &x_r, double &x_z)
Definition: basic_physics.hpp:139
KOKKOS_INLINE_FUNCTION void get_ExB_drift(double D, double drift_on, const SimdVector &bfield, const SimdVector &E, int i_simd, double over_B2, double(&yp_exb)[3])
Definition: push_physics.tpp:197
Simd< double > phi
Definition: simd.hpp:152
Simd< double > z
Definition: simd.hpp:151
Definition: simd.hpp:139
KOKKOS_INLINE_FUNCTION void get_curl_nb(const SimdVector &bfield, const double(&grad_B)[3], double over_B, double over_B2, const double(&curl_B)[3], double inv_r, double(&curl_nb)[3], int i_simd)
Definition: push_physics.tpp:141
phi coordinate
Definition: globals.hpp:203
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< GFPT > &fld) const
Definition: grid_field_pack.hpp:74
KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector &vec, int i_simd, const SimdVector &B, double Bmag) const
Definition: grid_field_pack.hpp:60
Definition: species.hpp:75
KOKKOS_INLINE_FUNCTION double magnitude(const int i_simd) const
Definition: simd.hpp:176
KOKKOS_INLINE_FUNCTION void get_v_par_drift(double D, const SimdVector &bfield, const SimdVector &tdb, const double(&curl_nb)[3], const LocalFields< GFPT > &fld, double over_B, double cmrho, double gradpsi_r_norm, double gradpsi_z_norm, double(&v_par)[3], int i_simd)
Definition: get_drift_velocity.hpp:62
z coordinate
Definition: globals.hpp:202
KOKKOS_INLINE_FUNCTION void fields_at_node(const GridFieldPack< DeviceType, GFPT > &gfpack, const PushControls &push_controls, const Grid< DeviceType > &grid, const SimdVector &B, const SimdVector2D &gradpsi, const Simd< int > &global_inds, const SimdGridWeights< Order::One, PIT_GLOBAL > grid_wts, const SimdGyroWeights< KinType::GyroKin > &rho_wts, LocalFields< GFPT > &fld)
Definition: get_drift_velocity.hpp:8