basic_physics.hpp

Go to the documentation of this file.

#ifndef BASIC_PHYSICS_HPP
#define BASIC_PHYSICS_HPP
#include "space_settings.hpp"
#include "constants.hpp"
#include "simd.hpp"

// Kinetic energy
// E = (1/2) m v^2
KOKKOS_INLINE_FUNCTION double kinetic_energy(double mass, double v){
    return 0.5*mass*v*v;
}

// Return equilibrium distribution function
// f = (n / sqrt(T^3)) * exp(-E/T)
// T and E are EV units for XGC normalization-particle.
KOKKOS_INLINE_FUNCTION double maxwellian_dist(double den, double temp, double energy){
    return den/(temp*sqrt(temp))*exp(-min(EXP_LIM,energy/temp));
}

// Maxwellian distribution of 1D velocity part
//  sqrt(m/TWOPI e) is missing when temp and energy are in eV unit.
//  sqrt(m/TWOPI)   is missing when temp and energy are in Joule.
// (1.0 - exp(-maxe)) is correction factor since the energy is finite
KOKKOS_INLINE_FUNCTION double maxwellian_dist_1d_v(double temp, double energy, double maxe){
    return 1.0/sqrt(temp)*exp(-min(EXP_LIM,energy/temp)) / (1.0-exp(-maxe));
}

// Maxwellian distribution of 2D velocity part
// m/(TWOPI e) is missing when temp and energy are in eV unit.
// m/TWOPI     is missing when temp and energy are in Joule.
// (1.0 - exp(-maxe)) is correction factor since the energy is finite
KOKKOS_INLINE_FUNCTION double maxwellian_dist_2d_v(double temp, double energy, double maxe){
    return 1.0/temp*exp(-min(EXP_LIM,energy/temp)) / (1.0-exp(-maxe));
}

// Return adiabatic exponential factor
// Should be the same as get_ad of f0module.F90
//    When the modified adiabatic function is used, the adiabatic response becomes
//    linear function in (-a,a) --> exp_ad(x) = 1+x
//    exponential decay to zero in x<-a --> exp_ad(x)=exp(-|x|+a)*(1-a)
//    This satisfies continuous condition
//    mirror symmetry (0,1) --> x-> -x, y-> 2-y
//     for (x>a), exp_ad(x)=2 - exp(-|x|+a)*(1-a)
KOKKOS_INLINE_FUNCTION double exp_ad(double x){
#ifdef USE_LINEAR_ADIABATIC_RESPONSE
    constexpr bool use_linear = true;
#else
    constexpr bool use_linear = false;
#endif

#ifdef LINEAR_ADIABATIC_RESPONSE_BD
    constexpr double lim=LINEAR_ADIABATIC_RESPONSE_BD;
#else
    constexpr double lim=0.7;
#endif
    if constexpr(use_linear){
        double tmp = exp(-abs(x)+lim)*(1-lim);
        double tmp2 = (x<-lim ? tmp : 1.0+x);
        return (x>lim ? 2.0-tmp : tmp2);
    } else {
        return exp(x);
    }
}

// gyro_radius
KOKKOS_INLINE_FUNCTION double gyro_radius(double B, double mu, double c2_2m){
    return sqrt(mu/(B*c2_2m));
}

/* Returns thermal velocity
 * @param[in]   mass     species mass
 * @param[in]   temp_ev  reference temperature
 * return thermal velocity
 */
KOKKOS_INLINE_FUNCTION double thermal_velocity(double mass, double temp_ev){
    return sqrt((temp_ev*EV_2_J)/mass);
}

/* Returns velocity normalized by thermal velocity
 * @param[in]   mass     species mass
 * @param[in]   temp_ev  reference temperature
 * @param[in]   v        input velocity
 * return normalized velocity
 */
KOKKOS_INLINE_FUNCTION double normalize_to_vth(double mass, double temp_ev, double v){
    return v*sqrt(mass/(temp_ev*EV_2_J));
}

/* Returns normalized parallel velocity
 * @param[in]   c_m      species charge over mass
 * @param[in]   mass     species mass
 * @param[in]   B        magnetic field magnitude
 * @param[in]   temp_ev  reference temperature
 * @param[in]   rho      gyroradius
 * return normalized parallel velocity
 */
KOKKOS_INLINE_FUNCTION double normalized_v_para(double c_m, double mass, double B, double temp_ev, double rho){
    return c_m*rho*B*sqrt(mass/(temp_ev*EV_2_J));
}

/* Returns normalized sqrt(mu)
 * @param[in]   B        magnetic field magnitude
 * @param[in]   temp_ev  particle/reference temperature
 * @param[out]  mu       mu
 * return normalized sqrt(mu)
*/
KOKKOS_INLINE_FUNCTION double normalized_sqrt_mu(double B, double temp_ev, double mu){
    return sqrt(mu*2.0*B/(temp_ev*EV_2_J));
}

// Convert from {rho, mu} to {energy, pitch}
KOKKOS_INLINE_FUNCTION void rho_mu_to_en_pitch(double B, double c_m, double c2_2m, double mass, double rho, double mu, double& en, double& pitch){
    en = mu*B + c2_2m*(rho*rho*B*B);
    double v_para = c_m*rho*B;
    double v_perp2 = 2.0*mu*B/mass;
    pitch=v_para/sqrt(v_para*v_para + v_perp2);
}

// Convert from {energy, pitch} to {rho, mu}
KOKKOS_INLINE_FUNCTION void en_pitch_to_rho_mu(double B, double c_m, double mass, double en, double pitch, double& rho, double& mu){
    rho = pitch*sqrt(2.0*en/mass)/(B*c_m);
    mu = max(0.0,(1.0-pitch*pitch)*en/B);
}

KOKKOS_INLINE_FUNCTION void rz_2_psitheta(double dpsidr, double dpsidz, double x_r, double x_z, double& x_rad, double& x_theta){
    x_rad   =  dpsidr * x_r + dpsidz * x_z;
    x_theta = -dpsidz * x_r + dpsidr * x_z;
}

KOKKOS_INLINE_FUNCTION void rz_2_psitheta(const Simd<double>& dpsidr, const Simd<double>& dpsidz, SimdVector2D& x){
    for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
        double x_rad, x_theta;
        rz_2_psitheta(dpsidr[i_simd], dpsidz[i_simd], x.r[i_simd], x.z[i_simd], x_rad, x_theta);
        x.r[i_simd] = x_rad;
        x.z[i_simd] = x_theta;
    }
}

KOKKOS_INLINE_FUNCTION void psitheta_2_rz(double dpsidr, double dpsidz, double x_rad, double x_theta, double& x_r, double& x_z){
    x_r = dpsidr * x_rad - dpsidz * x_theta;
    x_z = dpsidz * x_rad + dpsidr * x_theta;
}

#endif