CollisionModel.cpp

Go to the documentation of this file.

#include "CollisionModel.h"
#include "numeric/TableLookup.h"
#include "numeric/TrialStateError.h"
#include "units.h"
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <string>
 
namespace rizer {
 
namespace {
 
// Physical constants from units.h (mirror rizer/misc/units.py) so the C++
// collision physics reproduces the Python MixtureCollisionFrequencies
// 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;
 
// ---------------------------------------------------------------------------------- //
//   File-local helpers                                                               //
// ---------------------------------------------------------------------------------- //
 
// 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.transport.coulomb.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, refused where it is undefined exactly as the
// Python reference refuses it: without electrons the Debye length is
// infinite, and a Debye length below the mean impact parameter makes
// lnLambda negative, i.e. a negative collision frequency. A TrialStateError
// rather than a std:: exception: these depend on the state CVODE proposed,
// so it should retry a smaller step, not abandon the solve.
// @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.transport.coulomb.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)) {
        throw TrialStateError(
            "CollisionModel::coulombLogMitchner",
            "the Coulomb logarithm needs a positive electron "
            "density and temperature, got n_e = " + std::to_string(n_e)
            + " m^-3, Te = " + std::to_string(Te) + " K.");
    }
    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);
    if (lambda_D < b0_bar) {
        throw TrialStateError(
            "CollisionModel::coulombLogMitchner",
            "the Debye length " + std::to_string(lambda_D)
            + " m is below the mean impact parameter " + std::to_string(b0_bar)
            + " m (n_e = " + std::to_string(n_e) + " m^-3, Te = "
            + std::to_string(Te) + " K), so the Coulomb logarithm is undefined.");
    }
    return std::log(lambda_D / b0_bar);
}
 
// 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
 
// ---------------------------------------------------------------------------------- //
//   Constructor                                                                      //
// ---------------------------------------------------------------------------------- //
 
CollisionModel::CollisionModel(const std::vector<SpeciesSpec>& specs, double Te_min,
                               double Te_max, std::size_t Te_n)
{
    // 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) {
        throw std::invalid_argument(
            "CollisionModel: Te_n must be at least 2 to span [Te_min, Te_max], got "
            + std::to_string(Te_n) + ".");
    }
    if (!(Te_min > 0.0 && Te_max > Te_min)) {
        throw std::invalid_argument(
            "CollisionModel: the Te grid needs 0 < Te_min < Te_max, got Te_min = "
            + std::to_string(Te_min) + " K, Te_max = " + std::to_string(Te_max) + " K.");
    }
    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];
        SpeciesCollisionConfig& 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 = SpeciesCollisionKind::Skip;
        } else if (s.Z > 0) {
            sp.kind = SpeciesCollisionKind::Ion;
            sp.Z = s.Z;
            // Precomputed once here exactly as
            // IonMomentumTransferCollisionFrequencyModel.__init__ does in the
            // Python reference (rizer/transport/collision_frequency.py).
            sp.cst = ionCollisionPrefactor(s.Z);
        } else if (!s.energy_J.empty() || !s.sigma_m2.empty()) {
            // Validated, not used as a branch test. Letting a mismatched pair
            // fall through to the next arm silently swaps this species'
            // collision model for a hard sphere, or drops it from the heavy
            // sums entirely -- so a requested cross section would never be
            // used and nothing would say so.
            if (s.energy_J.size() != s.sigma_m2.size() || s.energy_J.size() < 2) {
                throw std::invalid_argument(
                    "CollisionModel: species " + std::to_string(k)
                    + " has a tabulated cross section of "
                    + std::to_string(s.energy_J.size()) + " energies against "
                    + std::to_string(s.sigma_m2.size())
                    + " values. They must be the same length and at least 2 "
                      "points. Leave both empty to use the hard-sphere radius "
                      "or to skip this species.");
            }
            sp.kind = SpeciesCollisionKind::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 = SpeciesCollisionKind::HardSphere;
            sp.Qbar = PropertyTable::constant(4.0 / 3.0 * kPI * s.radius * s.radius);
        } else {
            sp.kind = SpeciesCollisionKind::Skip;  // no model -> ignored in the heavy sums
        }
    }
}
 
// ---------------------------------------------------------------------------------- //
//   Public methods                                                                   //
// ---------------------------------------------------------------------------------- //
 
void CollisionModel::evaluate(double Te, double n_e, const double* n_k, std::size_t nsp,
                              double& sigma, double& nu_E, double& nu_eH, double& nu_eI,
                              double& nu_eH_mass_weighted) const
{
    sigma = 0.0;
    nu_E = 0.0;
    nu_eH = 0.0;
    nu_eI = 0.0;
    nu_eH_mass_weighted = 0.0;
    // Refused rather than answered with zeros: every frequency below carries a
    // power of Te, and sigma = nu_E = 0 would read as a cold neutral gas
    // instead of as a state with no electron temperature. The negated
    // comparison also rejects a NaN. A TrialStateError, so CVODE retries.
    if (!(Te > 0.0)) {
        throw TrialStateError(
            "CollisionModel::evaluate", "Te = " + std::to_string(Te)
            + " K; the electron temperature must be positive.");
    }
    if (m_species.empty()) return;
    // Refused rather than evaluated over the shorter of the two: that drops
    // species from every collision sum, or ignores the caller's extra ones.
    if (nsp != m_species.size()) {
        throw std::invalid_argument(
            "CollisionModel::evaluate: got number densities for "
            + std::to_string(nsp) + " species, but the model was built for "
            + std::to_string(m_species.size()) + ".");
    }
 
    const double v_th = thermalVelocity(Te);
    double nu_eH_sum = 0.0, nu_ions_sum = 0.0, nu_E_sum = 0.0, nu_eH_mass_weighted_sum = 0.0;
 
    for (std::size_t k = 0; k < nsp; k++) {
        const SpeciesCollisionConfig& sp = m_species[k];
        if (sp.kind == SpeciesCollisionKind::Skip) continue;
        const double n = n_k[k];
        if (n <= 0.0) continue;
 
        double nu1 = 0.0;
        if (sp.kind == SpeciesCollisionKind::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_sum += 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]
            if (!sp.Qbar.covers(Te)) {
                throwOutsideTable(sp.Qbar, Te, "CollisionModel::evaluate",
                                  "the collision cross section Qbar(Te) of species " +
                                      std::to_string(k),
                                  "K", "Widen the collision-model grid with Te_min/Te_max.");
            }
            nu1 = n * sp.Qbar.eval(Te) * v_th;
        }
        nu_eH_sum += nu1;
        // Elastic energy-exchange frequency contribution (Mitchner II-7.6),
        // the C++ counterpart of the same formula computed in Python by
        // rizer.transport.mixture_law.MixtureCollisionFrequencies.plasma_power_elastic:
        // @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
        // MixtureCollisionFrequencies.mass_weighted_electron_heavy_collision_frequency:
        // @f$\sum_h \bar\nu_{eh}/m_h@f$.
        nu_eH_mass_weighted_sum += sp.inv_mh * nu1;
    }
 
    nu_E = nu_E_sum;
    nu_eH = nu_eH_sum;
    nu_eI = nu_ions_sum;
    nu_eH_mass_weighted = nu_eH_mass_weighted_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.
        const double f_ion = nu_eI / nu_eH;
        sigma *= 1.0 + 0.98 * f_ion;
    }
}
 
// ---------------------------------------------------------------------------------- //
//   Public static methods                                                            //
// ---------------------------------------------------------------------------------- //
 
double CollisionModel::electronElectronCollisionFrequency(double n_e, double Te)
{
    // Mitchner II-8.11e with Z=1, the C++ counterpart of
    // MixtureCollisionFrequencies.electron_electron_collision_frequency: same
    // closed form as the ion collision frequency (ionCollisionPrefactor
    // above), evaluated with n=n_e, Z=1. No electrons is a real state with no
    // electron-electron collisions; a non-positive Te is not a state at all.
    if (!(Te > 0.0)) {
        throw TrialStateError(
            "CollisionModel::electronElectronCollisionFrequency",
            "Te = " + std::to_string(Te) + " K; the electron temperature must be positive.");
    }
    if (n_e <= 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;
}
 
 
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) {
        throw std::invalid_argument(
            "CollisionModel::meanCrossSection: the trapezoid rule needs N >= 2 "
            "points, got N = " + std::to_string(N) + ".");
    }
    if (energy_J.empty()) {
        throw std::invalid_argument(
            "CollisionModel::meanCrossSection: the cross section table is empty.");
    }
    // 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,
                                         "the cross section 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;
}
 
} // namespace rizer