get_drift_velocity.hpp

Go to the documentation of this file.

#ifndef GET_DRIFT_VELOCITY_HPP
#define GET_DRIFT_VELOCITY_HPP

#include "push_physics.hpp"

// Modified from fields_at_point in grid_field_pack. Should consolidate and use GridWeights<Order::Zero>
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){
    // Initialize local fields to 0
    fld.template get<Label::E>().zero();
#ifdef EXPLICIT_EM
    fld.template get<Label::dAh>().zero();
    fld.template get<Label::dAs>().zero();
    fld.template get<Label::Ah>().zero();
    fld.template get<Label::As>().zero();
    fld.template get<Label::Epar_em>().zero();
#endif

    for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
        // Field basis correction
        FieldCorrection field_correction(gradpsi, i_simd);
      
        constexpr double wp = 1.0; // No triangle weights since this is on a vertex
    
        field_correction.set(i_simd, grid.uses_rz_basis(global_inds[i_simd]), gradpsi);

        gfpack.gather_all_fields(push_controls, i_simd, global_inds[i_simd], wp, field_correction, grid_wts, rho_wts, fld);
    
#ifdef XGC1 
        // Correct phi component of E-field
        double Bmag = B.magnitude(i_simd);

        // Get phi from para
        gfpack.phi_from_para(fld.template get<Label::E>(), i_simd, B, Bmag);
#ifdef EXPLICIT_EM
        gfpack.phi_from_para(fld.template get<Label::dAh>(), i_simd, B, Bmag);
        gfpack.phi_from_para(fld.template get<Label::dAs>(), i_simd, B, Bmag);
#endif
#endif
    }
}

// Modified from fields_at_point in grid_field_pack. Should consolidate and use GridWeights<Order::Zero>
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){
    // Initialize local fields to 0
    pot_out = 0.0;
#ifdef XGC1
    for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
        int irho = rho_wts.irho(i_simd);
        constexpr double wp = 1.0; // No triangle weights since this is on a vertex
        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));
    }
#endif
}

// Parallel drift
// v_parallel = D*(B_star*dH/dv_p *1/(m*B) + v_p*dB_RMP)
// But we move the m/c*u_para*curl(b) part to the magnetic drift term.
// The remaining terms form the actual parallel motion along 
// the perturbed fieldn --> slightly different B_star
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){
    // B_star = B+A_s*curl(b) + grad(A)xB
    double b_star[3];
    b_star[PIR] = bfield.r[i_simd];
    b_star[PIZ] = bfield.z[i_simd];
    b_star[PIP] = bfield.phi[i_simd];
#ifdef EXPLICIT_EM
    const auto& As = fld.template get<Label::As>();
    const auto& dAs = fld.template get<Label::dAs>();
    constexpr int NEG_DA = -1; // dAs needs sign reversal!
    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;
    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;
    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;
#else
    // Remove radial component of b_star
    double rad, theta;
    rz_2_psitheta(gradpsi_r_norm, gradpsi_z_norm, b_star[PIR], b_star[PIZ], rad, theta);
    psitheta_2_rz(gradpsi_r_norm, gradpsi_z_norm, 0.0, theta, b_star[PIR], b_star[PIZ]);

    // Get radial component of tdb and add to b_star
    double tdb_r_mod, tdb_z_mod;
    rz_2_psitheta(gradpsi_r_norm, gradpsi_z_norm, tdb.r[i_simd], tdb.z[i_simd], rad, theta);
    psitheta_2_rz(gradpsi_r_norm, gradpsi_z_norm, rad, 0.0, tdb_r_mod, tdb_z_mod);
    b_star[PIR] += tdb_r_mod;
    b_star[PIZ] += tdb_z_mod;
#endif

    // In the EM formalism, the parallel flow is v_para*B_star, i.e., he parallel flow
    // contribution to the radial flux can come from the vr and dvr --->
    // Subtract the curl(B) contribution, which is accounted for in v_mag
    double prefac = D*cmrho;
    v_par[PIR] = prefac*b_star[PIR];
    v_par[PIZ] = prefac*b_star[PIZ];
    v_par[PIP] = prefac*b_star[PIP];
}

// The magnetic cross-field drifts
// Remove the (numerical flux-surface average of the magnetic drifts.
// This approximation accepts a small error because the velocity dependence
// in "D" is not taken into accound in v_gradb_avg and v_curv_avg.
//
// v_drift = D*(FxB) = D*(Bxgrad(H))/B^2  !Use all terms except E in grad(H)
// f = f_mr + f_pert = f_mr + f_es + f_em
// f_mr --> Mirror force
// f_es --> Electric field
// f_em --> Electromagnetic forces
// relevant force  f = f_mr + f_em
//                   = -mu/c * grad(B) + (v_para-c/m*Ah)*grad(Ah)
//  -mu/c*[grad(B)xB]/B^2 --> Apply cancellation law
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){
    // Convert to r,z
    double v_curv_r_fsa, v_curv_z_fsa, v_grad_B_r_fsa, v_grad_B_z_fsa;
    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);
    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);

    double over_c = 1.0/charge;

    // Grad B drift is Bxgrad(|B|)
    double grad_B_drift[3];
    grad_B_drift[PIR] = bfield.phi[i_simd] * grad_B[PIZ]       - bfield.z[i_simd]   * grad_B[PIP]*inv_r;
    grad_B_drift[PIZ] = bfield.r[i_simd]   * grad_B[PIP]*inv_r - bfield.phi[i_simd] * grad_B[PIR];
    grad_B_drift[PIP] = bfield.z[i_simd]   * grad_B[PIR]       - bfield.r[i_simd]   * grad_B[PIZ];

#ifdef EXPLICIT_EM
    double over_B = 1.0/B;
    const auto& Ah = fld.template get<Label::Ah>();
    const auto& dAh = fld.template get<Label::dAh>();
    constexpr int NEG_DA = -1; // dAh needs sign reversal!
    double fr = -mu*over_c*grad_B[PIR] + NEG_DA*u*dAh.r[i_simd];
    double fz = -mu*over_c*grad_B[PIZ] + NEG_DA*u*dAh.z[i_simd];
    double fp = -mu*over_c*grad_B[PIP] + NEG_DA*u*dAh.phi[i_simd];
#ifdef EXPLICIT_EM_NONLIN_TERMS
    fr -= NEG_DA*c_m*Ah[i_simd]*dAh.r[i_simd];
    fz -= NEG_DA*c_m*Ah[i_simd]*dAh.z[i_simd];
    fp -= NEG_DA*c_m*Ah[i_simd]*dAh.phi[i_simd];
#endif

    // Magnetic drifts in (R,phi,Z) coordinates:
    v_mag[PIR] = D*(bfield.z[i_simd]*fp - bfield.phi[i_simd]*fz)*over_B2;
    v_mag[PIZ] = D*(bfield.phi[i_simd]*fr - bfield.r[i_simd]*fp)*over_B2;
    v_mag[PIP] = D*(bfield.r[i_simd]*fz - bfield.z[i_simd]*fr)*over_B2;

    // Apply cancellation law of Bxgrad(B)
    // Also add contribution from curl(b) (removed from B*dH/du)
    v_mag[PIR] -= D*over_c*over_B2*mu*v_grad_B_r_fsa;
    v_mag[PIZ] -= D*over_c*over_B2*mu*v_grad_B_z_fsa;

    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));
    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));
    v_mag[PIP] += D*over_B*u*rho*(over_B* grad_B_drift[PIP]                   +  B*curl_B[PIP]);
#else
    double cmrho2 = cmrho*rho;
    double murho2B = (mu + charge*cmrho2*B);
    double murho2B_c = murho2B*over_c;

    // The magnetic drifts (curvature + grad_b)
    // Remove the (numerical flux-surface average of the magnetic drifts.
    // This approximation accepts a small error because the velocity dependence
    // in "D" is not taken into accound in v_gradb_avg and v_curv_avg.
    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);
    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);
    v_mag[PIP] = D*( grad_B_drift[PIP]                  *murho2B_c*over_B2 +  curl_B[PIP]                *cmrho2);
#endif
}


//> Evaluates the drift velocities at given inode, imu, ivp
// The basic magnetic field quantities necessary for the calculation of the
// magnetic drifts are stored in the grid data structure.
// To save memory, this could also be shifted to f0_module.
//
// CAUTION: Although this routine uses the full expression for dR/dt,
//          it is correct only when using the pullback method and when
//          in addition, A_h=0, i.e., directly after the pullback operation
//          because it uses the assumption that u_para=v_para
//
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){

    for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
        const double& B = B_mag[i_simd];
        double rho = vp[i_simd]/(B*species.c_m);

        double over_B = 1.0/B;
        double over_B2 = over_B*over_B;
        double cmrho = species.c_m*rho;

        // For clarity, distinguish v_para and u_para
        double u = vp[i_simd];
#ifdef EXPLICIT_EM
        u += species.c_m*fld.template get<Label::Ah>()[i_simd];
#endif

        double grad_B[3];
        get_grad_B(bfield, jacb, over_B, grad_B, i_simd);

        double curl_B[3];
        get_curl_B(bfield, jacb, inv_r[i_simd], curl_B, i_simd);
        
        double curl_nb[3]; // [ curl of normalized vec B ]  =  curl vec B / B + vec B X grad B / B^2
        get_curl_nb(bfield, grad_B, over_B, over_B2, curl_B, inv_r[i_simd], curl_nb, i_simd);

        // Velocity space Jacobian factor
#ifdef EXPLICIT_EM
        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);
#else       
        const double D = get_vspace_jacb_factor(rho, over_B2, bfield, jacb, inv_r[i_simd], i_simd);
#endif

        // The ExB drift
        constexpr double drift_on = 1.0;
        double yp[3];
        get_ExB_drift(D, drift_on, bfield, fld.template get<Label::E>(), i_simd, over_B2, inv_r[i_simd], yp);
        v_exb.r[i_simd] = yp[PIR];
        v_exb.z[i_simd] = yp[PIZ];
        v_exb.phi[i_simd] = yp[PIP]/inv_r[i_simd]; // Seems like a bug to have inv_r here? ALS

        // Magnetic drifts
        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);
        v_mag.r[i_simd] = yp[PIR];
        v_mag.z[i_simd] = yp[PIZ];
        v_mag.phi[i_simd] = yp[PIP];

        // Parallel drift
        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);
        v_par.r[i_simd] = yp[PIR];
        v_par.z[i_simd] = yp[PIZ];
        v_par.phi[i_simd] = yp[PIP];
    }
}

#endif