1 #ifndef GET_DRIFT_VELOCITY_HPP
2 #define GET_DRIFT_VELOCITY_HPP
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){
10 fld.template get<Label::E>().zero();
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();
19 for (
int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
23 constexpr
double wp = 1.0;
25 field_correction.
set(i_simd, grid.
uses_rz_basis(global_inds[i_simd]), gradpsi);
27 gfpack.
gather_all_fields(push_controls, i_simd, global_inds[i_simd], wp, field_correction, grid_wts, rho_wts, fld);
34 gfpack.
phi_from_para(fld.template get<Label::E>(), i_simd, B, Bmag);
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);
44 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){
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;
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));
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){
65 b_star[
PIR] = bfield.
r[i_simd];
66 b_star[
PIZ] = bfield.
z[i_simd];
67 b_star[
PIP] = bfield.
phi[i_simd];
69 const auto&
As = fld.template get<Label::As>();
70 const auto&
dAs = fld.template get<Label::dAs>();
71 constexpr
int NEG_DA = -1;
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;
76 b_star[
PIR] += tdb.
r[i_simd];
77 b_star[
PIZ] += tdb.
z[i_simd];
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];
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){
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);
109 double over_c = 1.0/charge;
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];
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;
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];
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;
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;
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]);
145 double cmrho2 = cmrho*rho;
146 double murho2B = (mu + charge*cmrho2*B);
147 double murho2B_c = murho2B*over_c;
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);
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){
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);
177 double over_B = 1.0/B;
178 double over_B2 = over_B*over_B;
179 double cmrho = species.
c_m*rho;
182 double u = vp[i_simd];
184 u += species.
c_m*fld.template get<Label::Ah>()[i_simd];
188 get_grad_B(bfield, jacb, over_B, grad_B, i_simd);
191 get_curl_B(bfield, jacb, inv_r[i_simd], curl_B, i_simd);
194 get_curl_nb(bfield, grad_B, over_B, over_B2, curl_B, inv_r[i_simd], curl_nb, i_simd);
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);
204 constexpr
double drift_on = 1.0;
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];
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];
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];
KOKKOS_INLINE_FUNCTION void set(int i_simd, double basis, const SimdVector2D &gradpsi)
Definition: field.hpp:35
Simd< double > r
Definition: simd.hpp:150
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
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
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