CollisionModel.cpp

Go to the documentation of this file.

#include "CollisionModel.h"
#include "units.h"
 
#include <algorithm>
#include <cmath>
 
namespace rizer {
 
namespace {
 
// Physical constants from units.h (mirror rizer/misc/units.py) so the C++
// collision physics reproduces the Python PlasmaExtension reference to
// interpolation error. Short local aliases keep the formula bodies readable.
constexpr double kE   = units::e;          // elementary charge [C]
constexpr double kKB  = units::k_b;        // Boltzmann constant [J/K]
constexpr double kME   = units::m_e;       // electron mass [kg]
constexpr double kEPS0 = units::epsilon_0; // vacuum permittivity [F/m]
constexpr double kNA   = units::N_a;       // Avogadro [1/mol]
constexpr double kPI   = units::pi;
 
// Electron thermal velocity (equations.py):
// @f[
//     v_\text{th} = \sqrt{\frac{8 k_B T_e}{\pi m_e}} \quad [\text{m/s}].
// @f]
// C++ counterpart of rizer.plasma.equations.electron_thermal_velocity.
// \param Te  Electron temperature [K]
inline double thermalVelocity(double Te)
{
    return std::sqrt(8.0 * kKB * Te / (kPI * kME));
}
 
// Mitchner Coulomb logarithm, floored at 1 for robustness during violent
// transients (the Python model raises for lnLambda <= 0 and warns for <= 1;
// a small positive floor keeps nu_ei finite and positive so the stiff
// solver never sees a non-physical collision frequency):
// @f[
//     \ln\Lambda = \ln\!\left(\frac{\lambda_D}{\overline{b_0}}\right), \qquad
//     \lambda_D = \sqrt{\frac{\varepsilon_0 k_B T_e}{n_e e^2}}, \qquad
//     \overline{b_0} = \frac{Z e^2}{12 \pi \varepsilon_0 k_B T_e}.
// @f]
// C++ counterpart of rizer.plasma.equations.coulomb_logarithm(model="Mitchner").
// \param n_e  Electron number density [m^-3]
// \param Te  Electron temperature [K]
// \param Z  Ion charge state (>0) of the colliding heavy species
inline double coulombLogMitchner(double n_e, double Te, int Z)
{
    if (n_e <= 0.0 || Te <= 0.0) return 1.0;
    const double lambda_D = std::sqrt(kEPS0 * kKB * Te / (n_e * kE * kE));
    const double b0_bar = static_cast<double>(Z) * kE * kE
                          / (12.0 * kPI * kEPS0 * kKB * Te);
    const double lnL = std::log(lambda_D / b0_bar);
    return (lnL > 1.0) ? lnL : 1.0;
}
 
// Mitchner II-8.11e prefactor (includes Z^2):
// @f[
//     c_\text{st}(Z) = \frac{4 \sqrt{2 \pi}}{3}
//         \left(\frac{m_e}{k_B}\right)^{1.5}
//         \left(\frac{e^2}{4 \pi \varepsilon_0 m_e}\right)^2 Z^2.
// @f]
// Shared by the per-ion-species prefactor (CollisionModel ctor) and the
// electron-electron collision frequency (Z=1) below.
// \param Z  Charge state (ion charge number, or 1 for electron-electron)
double ionCollisionPrefactor(int Z)
{
    const double a = 4.0 * std::sqrt(2.0 * kPI) / 3.0;
    const double b = std::pow(kME / kKB, 1.5);
    const double c = kE * kE / (4.0 * kPI * kEPS0 * kME);
    return a * b * c * c * static_cast<double>(Z) * static_cast<double>(Z);
}
 
} // namespace
 
double CollisionModel::meanCrossSection(const std::vector<double>& energy_J,
                                        const std::vector<double>& sigma_m2, double Te,
                                        double x_max, int N)
{
    // Trapezoid quadrature on a uniform x-grid of N points for (Mitchner
    // II-6.30):
    // @f[
    //     \overline{Q}(T_e) = \frac{2}{3} \int_0^{x_\text{max}}
    //         x^2 e^{-x} Q(k_B T_e x)\, dx,
    // @f]
    // identical to the Python TabulatedSpeciesCrossSection
    // .get_mean_cross_section (x_max=20, N=1000), which is the C++
    // counterpart of this method.
    if (N < 2 || energy_J.empty()) return 0.0;
    // np.interp default: linear interpolation, holding the nearest edge value
    // outside the table (constant extrapolation) rather than zeroing -- exactly
    // PropertyTable::eval()'s contract. Hard-zeroing biases the averaged cross
    // section -- and hence k(Te) -- low when Te samples energies beyond the
    // tabulated range. Built once here (not per query) since energy_J/sigma_m2
    // are fixed across this loop.
    const PropertyTable sigma_table = PropertyTable::fromTableOrScalar(energy_J, sigma_m2);
    const double dx = x_max / static_cast<double>(N - 1);
    double integral = 0.0;
    double prev = 0.0;  // integrand at x=0 is 0 (x^2 factor)
    for (int i = 1; i < N; i++) {
        const double x = dx * static_cast<double>(i);
        const double energy = kKB * Te * x;
        const double Q = sigma_table.eval(energy);
        const double cur = (2.0 / 3.0) * x * x * std::exp(-x) * Q;
        integral += 0.5 * (prev + cur) * dx;
        prev = cur;
    }
    return integral;
}
 
std::vector<double> CollisionModel::meanCrossSectionTable(
    const std::vector<double>& energy_J, const std::vector<double>& sigma_m2,
    const std::vector<double>& Te, double x_max, int N)
{
    std::vector<double> out(Te.size());
    for (std::size_t i = 0; i < Te.size(); i++) {
        out[i] = meanCrossSection(energy_J, sigma_m2, Te[i], x_max, N);
    }
    return out;
}
 
CollisionModel::CollisionModel(const std::vector<SpeciesSpec>& specs, double Te_min,
                               double Te_max, std::size_t Te_n, bool spitzer_correction)
    : m_spitzer(spitzer_correction)
{
    // Log-spaced Te grid for the precomputed Qbar(Te) tables (smooth function;
    // log spacing resolves the low-Te knee). PropertyTable handles non-uniform.
    if (Te_n < 2) Te_n = 2;
    std::vector<double> Te_grid(Te_n);
    const double llo = std::log(Te_min), lhi = std::log(Te_max);
    for (std::size_t i = 0; i < Te_n; i++) {
        Te_grid[i] = std::exp(llo + (lhi - llo) * static_cast<double>(i)
                                            / static_cast<double>(Te_n - 1));
    }
 
    m_species.resize(specs.size());
    for (std::size_t k = 0; k < specs.size(); k++) {
        const SpeciesSpec& s = specs[k];
        Sp& sp = m_species[k];
        if (s.molar_mass > 0.0) {
            // m_h [kg] = molar_mass [kg/mol] / N_A, so N_A / molar_mass = 1 / m_h;
            // molar_mass is already SI (kg/mol) here -- the Python side converts
            // from g/mol via *1e-3 before this ever reaches C++ (see
            // unpack_momentum_transfer_models / PlasmaChannel species setup).
            sp.two_me_over_mh = 2.0 * kME * kNA / s.molar_mass;  // 2 m_e / m_h
            sp.inv_mh = kNA / s.molar_mass;                      // 1 / m_h
        }
        if (s.is_electron) {
            sp.kind = Kind::Skip;
        } else if (s.Z > 0) {
            sp.kind = Kind::Ion;
            sp.Z = s.Z;
            // Precomputed once here exactly as
            // IonMomentumTransferCollisionFrequencyModel.__init__ does in the
            // Python reference (rizer/plasma/collision_frequency.py).
            sp.cst = ionCollisionPrefactor(s.Z);
        } else if (!s.energy_J.empty() && s.energy_J.size() == s.sigma_m2.size()) {
            sp.kind = Kind::Tabulated;
            // C++ counterpart of TabulatedSpeciesCrossSection: bake the Maxwellian
            // average into a Te table once at construction (build time), rather
            // than integrating live every RHS evaluation.
            std::vector<double> Qbar =
                meanCrossSectionTable(s.energy_J, s.sigma_m2, Te_grid);
            sp.Qbar = PropertyTable(Te_grid, Qbar);
        } else if (s.radius > 0.0) {
            sp.kind = Kind::HardSphere;
            sp.Qbar = PropertyTable::constant(4.0 / 3.0 * kPI * s.radius * s.radius);
        } else {
            sp.kind = Kind::Skip;  // no model -> ignored in the heavy sums
        }
    }
}
 
void CollisionModel::evaluate(double Te, double n_e, const double* n_k, std::size_t nsp,
                              double& sigma, double& nu_E, double& cond_maxwell_2) const
{
    sigma = 0.0;
    nu_E = 0.0;
    cond_maxwell_2 = 0.0;
    if (m_species.empty() || Te <= 0.0) return;
 
    const double v_th = thermalVelocity(Te);
    double nu_eH = 0.0, nu_ions = 0.0, nu_E_sum = 0.0, cond_maxwell_2_sum = 0.0;
 
    const std::size_t K = std::min(nsp, m_species.size());
    for (std::size_t k = 0; k < K; k++) {
        const Sp& sp = m_species[k];
        if (sp.kind == Kind::Skip) continue;
        const double n = n_k[k];
        if (n <= 0.0) continue;
 
        double nu1 = 0.0;
        if (sp.kind == Kind::Ion) {
            // Mitchner II-8.10, the C++ counterpart of
            // IonMomentumTransferCollisionFrequencyModel
            // .get_mean_momentum_transfer_collision_frequency:
            // @f[
            //     \nu_{ei} = n_i \, c_\text{st} \, T_e^{-1.5} \, \ln\Lambda.
            // @f]
            const double lnL = coulombLogMitchner(n_e, Te, sp.Z);
            nu1 = n * sp.cst * std::pow(Te, -1.5) * lnL;
            nu_ions += nu1;
        } else {  // Tabulated or HardSphere
            // Mitchner II-6.29, the C++ counterpart of
            // MomentumTransferCollisionFrequencyModel
            // .get_mean_momentum_transfer_collision_frequency:
            // @f[
            //     \nu_{eh} = n_h \, \overline{Q}(T_e) \, v_\text{th}.
            // @f]
            nu1 = n * sp.Qbar.eval(Te) * v_th;
        }
        nu_eH += nu1;
        // Elastic energy-exchange frequency contribution (Mitchner II-7.6),
        // the C++ counterpart of EnergyExchangeCollisionFrequencyModel
        // .get_mean_energy_exchange_collision_frequency:
        // @f[
        //     \nu_{eh}^E = \frac{2 m_e}{m_h} \nu_{eh}.
        // @f]
        nu_E_sum += sp.two_me_over_mh * nu1;
        // Second Maxwellian-distribution condition term (Mitchner VIII-3.8),
        // the C++ counterpart of
        // PlasmaExtension.conditions_for_maxwellian_distribution's
        // cond_maxwell_2: @f$\sum_h \bar\nu_{eh}/m_h@f$.
        cond_maxwell_2_sum += sp.inv_mh * nu1;
    }
 
    nu_E = nu_E_sum;
    cond_maxwell_2 = cond_maxwell_2_sum;
    if (nu_eH > 0.0 && n_e > 0.0) {
        sigma = n_e * kE * kE / (kME * nu_eH);
        // Electron-electron collisions raise sigma above the Lorentz value by up to
        // the Spitzer factor 1.98 for a fully ionized, singly charged plasma
        // (Mitchner II-13.18). Bridge it smoothly with the ion collision fraction
        // @f$f_\text{ion} = \nu_\text{ions} / \nu_{eH}@f$ (1 -> neutral-dominated,
        // 1.98 -> fully ionized) instead of switching abruptly once nu_ions
        // exceeds nu_neutrals.
        if (m_spitzer) {
            const double f_ion = nu_ions / nu_eH;
            sigma *= 1.0 + 0.98 * f_ion;
        }
    }
}
 
double CollisionModel::conductivity(double Te, double n_e, const double* n_k,
                                    std::size_t nsp) const
{
    double sigma = 0.0, nu_E = 0.0, cond_maxwell_2 = 0.0;
    evaluate(Te, n_e, n_k, nsp, sigma, nu_E, cond_maxwell_2);
    return sigma;
}
 
double CollisionModel::electronElectronCollisionFrequency(double n_e, double Te)
{
    // Mitchner II-8.11e with Z=1, the C++ counterpart of
    // PlasmaExtension.conditions_for_maxwellian_distribution's nu_ee: same
    // closed form as the ion collision frequency (ionCollisionPrefactor
    // above), evaluated with n=n_e, Z=1.
    if (n_e <= 0.0 || Te <= 0.0) return 0.0;
    static const double cst_ee = ionCollisionPrefactor(1);
    const double lnL = coulombLogMitchner(n_e, Te, 1);
    return n_e * cst_ee * std::pow(Te, -1.5) * lnL;
}
 
} // namespace rizer