1 #ifndef GET_DRIFT_VELOCITY_HPP
2 #define GET_DRIFT_VELOCITY_HPP
7 KOKKOS_INLINE_FUNCTION
void fields_at_node(
const GridFieldPack<DeviceType, PIT_GLOBAL>& 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& fld){
9 fld.template get<Label::E>().zero();
11 fld.template get<Label::dAh>().zero();
12 fld.template get<Label::dAs>().zero();
13 fld.template get<Label::Ah>().zero();
14 fld.template get<Label::As>().zero();
15 fld.template get<Label::Epar_em>().zero();
18 for (
int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
22 constexpr
double wp = 1.0;
24 field_correction.
set(i_simd, grid.
uses_rz_basis(global_inds[i_simd]), gradpsi);
26 gfpack.
gather_all_fields(push_controls, i_simd, global_inds[i_simd], wp, field_correction, grid_wts, rho_wts, fld);
33 gfpack.
phi_from_para(fld.template get<Label::E>(), i_simd, B, Bmag);
35 gfpack.
phi_from_para(fld.template get<Label::dAh>(), i_simd, B, Bmag);
36 gfpack.
phi_from_para(fld.template get<Label::dAs>(), i_simd, B, Bmag);
43 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){
47 for (
int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
48 int irho = rho_wts.irho(i_simd);
49 constexpr
double wp = 1.0;
50 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));
60 KOKKOS_INLINE_FUNCTION
void get_v_par_drift(
double D,
const SimdVector& bfield,
const SimdVector& tdb,
const double (&curl_nb)[3],
const LocalFields& fld,
double over_B,
double cmrho,
double gradpsi_r_norm,
double gradpsi_z_norm,
double (&v_par)[3],
int i_simd){
63 b_star[
PIR] = bfield.
r[i_simd];
64 b_star[
PIZ] = bfield.
z[i_simd];
65 b_star[
PIP] = bfield.
phi[i_simd];
67 const auto&
As = fld.template get<Label::As>();
68 const auto&
dAs = fld.template get<Label::dAs>();
69 constexpr
int NEG_DA = -1;
70 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;
71 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;
72 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;
74 b_star[
PIR] += tdb.
r[i_simd];
75 b_star[
PIZ] += tdb.
z[i_simd];
81 double prefac = D*cmrho;
82 v_par[
PIR] = prefac*b_star[
PIR];
83 v_par[
PIZ] = prefac*b_star[
PIZ];
84 v_par[
PIP] = prefac*b_star[
PIP];
100 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& 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){
102 double v_curv_r_fsa, v_curv_z_fsa, v_grad_B_r_fsa, v_grad_B_z_fsa;
103 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);
104 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);
106 double over_c = 1.0/charge;
109 double grad_B_drift[3];
110 grad_B_drift[
PIR] = bfield.
phi[i_simd] * grad_B[
PIZ] - bfield.
z[i_simd] * grad_B[
PIP]*inv_r;
111 grad_B_drift[
PIZ] = bfield.
r[i_simd] * grad_B[
PIP]*inv_r - bfield.
phi[i_simd] * grad_B[
PIR];
112 grad_B_drift[
PIP] = bfield.
z[i_simd] * grad_B[
PIR] - bfield.
r[i_simd] * grad_B[
PIZ];
115 double over_B = 1.0/B;
116 const auto&
Ah = fld.template get<Label::Ah>();
117 const auto&
dAh = fld.template get<Label::dAh>();
118 constexpr
int NEG_DA = -1;
119 double fr = -mu*over_c*grad_B[
PIR] + NEG_DA*u*dAh.r[i_simd];
120 double fz = -mu*over_c*grad_B[
PIZ] + NEG_DA*u*dAh.z[i_simd];
121 double fp = -mu*over_c*grad_B[
PIP] + NEG_DA*u*dAh.phi[i_simd];
122 #ifdef EXPLICIT_EM_NONLIN_TERMS
123 fr -= NEG_DA*c_m*
Ah[i_simd]*dAh.r[i_simd];
124 fz -= NEG_DA*c_m*
Ah[i_simd]*dAh.z[i_simd];
125 fp -= NEG_DA*c_m*
Ah[i_simd]*dAh.phi[i_simd];
129 v_mag[
PIR] = D*(bfield.
z[i_simd]*fp - bfield.
phi[i_simd]*fz)*over_B2;
130 v_mag[
PIZ] = D*(bfield.
phi[i_simd]*fr - bfield.
r[i_simd]*fp)*over_B2;
131 v_mag[
PIP] = D*(bfield.
r[i_simd]*fz - bfield.
z[i_simd]*fr)*over_B2;
135 v_mag[
PIR] -= D*over_c*over_B2*mu*v_grad_B_r_fsa;
136 v_mag[
PIZ] -= D*over_c*over_B2*mu*v_grad_B_z_fsa;
138 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));
139 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));
140 v_mag[
PIP] += D*over_B*u*rho*(over_B* grad_B_drift[
PIP] + B*curl_B[
PIP]);
142 double cmrho2 = cmrho*rho;
143 double murho2B = (mu + charge*cmrho2*B);
144 double murho2B_c = murho2B*over_c;
150 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);
151 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);
152 v_mag[
PIP] = D*( grad_B_drift[
PIP] *murho2B_c*over_B2 + curl_B[
PIP] *cmrho2);
167 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& fld,
SimdVector& v_mag,
SimdVector& v_exb,
SimdVector& v_par){
169 for (
int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
170 const double& B = B_mag[i_simd];
171 double rho = vp[i_simd]/(B*species.
c_m);
173 double over_B = 1.0/B;
174 double over_B2 = over_B*over_B;
175 double cmrho = species.
c_m*rho;
178 double u = vp[i_simd];
180 u += species.
c_m*fld.template get<Label::Ah>()[i_simd];
184 get_grad_B(bfield, jacb, over_B, grad_B, i_simd);
187 get_curl_B(bfield, jacb, inv_r[i_simd], curl_B, i_simd);
190 get_curl_nb(bfield, grad_B, over_B, over_B2, curl_B, inv_r[i_simd], curl_nb, i_simd);
194 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);
200 constexpr
double drift_on = 1.0;
202 get_ExB_drift(D, drift_on, bfield, fld.template get<Label::E>(), i_simd, over_B2, yp);
203 v_exb.
r[i_simd] = yp[
PIR];
204 v_exb.
z[i_simd] = yp[
PIZ];
205 v_exb.
phi[i_simd] = yp[
PIP];
208 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);
209 v_mag.
r[i_simd] = yp[
PIR];
210 v_mag.
z[i_simd] = yp[
PIZ];
211 v_mag.
phi[i_simd] = yp[
PIP];
214 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);
215 v_par.
r[i_simd] = yp[
PIR];
216 v_par.
z[i_simd] = yp[
PIZ];
217 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
double c_m
c/m
Definition: species.hpp:86
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 &fld) const
Definition: grid_field_pack.hpp:127
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:21
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:251
Definition: local_fields.hpp:33
Definition: grid_field.hpp:22
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 &fld, SimdVector &v_mag, SimdVector &v_exb, SimdVector &v_par)
Definition: get_drift_velocity.hpp:167
Definition: gyro_radius.hpp:24
r coordinate
Definition: globals.hpp:186
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 &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:100
KOKKOS_INLINE_FUNCTION void phi_from_para(SimdVector &vec, int i_simd, const SimdVector &B, double Bmag) const
Definition: grid_field_pack.hpp:121
double charge
Particle charge.
Definition: species.hpp:84
KOKKOS_INLINE_FUNCTION void get_v_par_drift(double D, const SimdVector &bfield, const SimdVector &tdb, const double(&curl_nb)[3], const LocalFields &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:60
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:43
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:196
Simd< double > phi
Definition: simd.hpp:152
Simd< double > z
Definition: simd.hpp:151
KOKKOS_INLINE_FUNCTION void fields_at_node(const GridFieldPack< DeviceType, PIT_GLOBAL > &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 &fld)
Definition: get_drift_velocity.hpp:7
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:188
Definition: species.hpp:75
KOKKOS_INLINE_FUNCTION double magnitude(const int i_simd) const
Definition: simd.hpp:176
z coordinate
Definition: globals.hpp:187