PlasmaRates.cpp

Go to the documentation of this file.

// PlasmaRates.cpp — see PlasmaRates.h for the design rationale.
 
#include "PlasmaRates.h"
#include "units.h"
#include "numeric/TrialStateError.h"
#include "thermo/PlasmaState.h"
 
#include "cantera/kinetics/ReactionRateFactory.h"
#include "cantera/kinetics/Arrhenius.h"
#include "cantera/kinetics/MultiRate.h"
#include "cantera/kinetics/Reaction.h"
#include "cantera/kinetics/Kinetics.h"
#include "cantera/kinetics/KineticsFactory.h"
#include "cantera/thermo/PlasmaPhase.h"
#include "cantera/thermo/IdealGasPhase.h"
#include "cantera/base/Solution.h"
#include "cantera/base/ctexceptions.h"
 
#include <cmath>
#include <stdexcept>
#include <string>
 
using namespace Cantera;
 
namespace rizer {
 
// ---------------------------------------------------------------------------------- //
//   File-local helpers                                                               //
// ---------------------------------------------------------------------------------- //
 
namespace {

double logForwardRateConstant(const ReactionRate& rate, double kf, double T,
                              const std::string& equation)
{
    const std::string type = rate.type();
    if (type == "Arrhenius" || type == "two-temperature-plasma") {
        const auto& arrhenius = dynamic_cast<const ArrheniusBase&>(rate);
        const double A = arrhenius.preExponentialFactor();
        if (!(A > 0.0)) {
            throw std::domain_error(
                "logForwardRateConstant: forward reaction '" + equation +
                "' has a non-positive pre-exponential factor, which has no log.");
        }
        return std::log(A) + arrhenius.temperatureExponent() * std::log(T) -
               arrhenius.activationEnergy() / (GasConstant * T);
    }
    if (type == "Druyvesteyn") {
        return dynamic_cast<const DruyvesteynRate&>(rate).logRate(T);
    }
    if (kf > 0.0 && std::isfinite(kf)) {
        return std::log(kf);
    }
    throw std::domain_error(
        "logForwardRateConstant: the forward rate constant of '" + equation +
        "' (type " + type + ") is " + std::to_string(kf) + " at T = " +
        std::to_string(T) + " K, and that rate type has no log form to recover it.");
}

bool anyNonNormal(const std::vector<double>& kf, const std::vector<double>& K)
{
    for (std::size_t i = 0; i < kf.size(); i++) {
        if (!std::isnormal(kf[i]) || !std::isnormal(K[i])) {
            return true;
        }
    }
    return false;
}
 
} // namespace
 
// ---------------------------------------------------------------------------------- //
//   Registration                                                                     //
// ---------------------------------------------------------------------------------- //
 
std::shared_ptr<Cantera::Solution> loadMechanism(const std::string& mech,
                                                 const std::string& phase)
{
    registerPlasmaRates();  // Idempotent
    return Cantera::newSolution(mech, phase);
}
 
void registerPlasmaRates()
{
    auto* f = ReactionRateFactory::factory();
    // `reg` replaces any existing entry silently, so this is idempotent.
    f->reg("janev-dissociative-recombination-C2Hy",
        // Lambda parameters (a plain, non-Doxygen comment: this lambda has no
        // separate declaration for \param to attach to, so a doc-comment
        // here would be misattached to the enclosing registerPlasmaRates()):
        //   node   YAML/AnyMap node for this reaction's rate block
        //          (holds "A").
        //   units  Unit context used to convert "A" into Cantera's
        //          internal rate-coeff units.
        [](const AnyMap& node, const UnitStack& units)
        {
            return new JanevDissociativeRecombinationC2HyRate(node, units);
        });
    f->reg("Druyvesteyn",
        // Lambda parameters (plain comment, see the note above):
        //   node   YAML/AnyMap node for this reaction's rate block
        //          (holds "A", "b", "Ea").
        //   units  Unit context used to convert "A" into Cantera's
        //          internal rate-coeff units.
        [](const AnyMap& node, const UnitStack& units)
        {
            return new DruyvesteynRate(node, units);
        });
    f->reg("reverse-two-temperature-plasma",
        // Lambda parameters (plain comment, see the note above):
        //   node   YAML/AnyMap node for this reaction's rate block
        //          (holds "forward_equation", "T").
        //   units  Unit context passed through to
        //          ReactionRate::setParameters.
        [](const AnyMap& node, const UnitStack& units)
        {
            return new ReverseTwoTemperaturePlasmaRate(node, units);
        });
}
 
// ---------------------------------------------------------------------------------- //
//   janev-dissociative-recombination-C2Hy - shared data                              //
// ---------------------------------------------------------------------------------- //
 
bool ElectronTemperatureData::update(const ThermoPhase& phase, const Kinetics& /*kin*/)
{
    auto plasma = dynamic_cast<const PlasmaPhase*>(&phase);
    double Te = plasma ? plasma->electronTemperature() : phase.temperature();
    if (Te != electronTemp) {
        electronTemp = Te;
        return true;
    }
    return false;
}
 
// ---------------------------------------------------------------------------------- //
//   janev-dissociative-recombination-C2Hy - constructor                              //
// ---------------------------------------------------------------------------------- //
 
JanevDissociativeRecombinationC2HyRate::JanevDissociativeRecombinationC2HyRate(
    const Cantera::AnyMap& node, const Cantera::UnitStack& rate_units)
{
    setParameters(node, rate_units);
}
 
// ---------------------------------------------------------------------------------- //
//   janev-dissociative-recombination-C2Hy - public methods                           //
// ---------------------------------------------------------------------------------- //
 
std::unique_ptr<MultiRateBase>
JanevDissociativeRecombinationC2HyRate::newMultiRate() const
{
    return std::make_unique<
        MultiRate<JanevDissociativeRecombinationC2HyRate, ElectronTemperatureData>>();
}
 
void JanevDissociativeRecombinationC2HyRate::setParameters(
    const AnyMap& node, const UnitStack& rate_units)
{
    ReactionRate::setParameters(node, rate_units);  // sets units + m_input
    if (!node.hasKey("A")) {
        m_A = NAN;
        return;
    }
    // A has the same units as the rate constant (the Te-dependent factor is
    // dimensionless), so it converts like any rate coefficient (e.g.
    // cm^3/mol/s or m^3/kmol/s, per the reaction's declared `units:` block,
    // both into Cantera's internal kmol/m/s system).
    m_A = node.units().convertRateCoeff(node["A"], conversionUnits());
    m_valid = true;
}
 
// Compiled counterpart of the Python `JanevDissociativeRecombinationC2Hy`
// ct.ExtensibleRate (rizer/kin/extensible_rate.py), which implements
// Janev2004 Eq. 70 for dissociative recombination: a closed-form rate that
// depends on the electron temperature Te only (no Tg, no reduced field).
double JanevDissociativeRecombinationC2HyRate::evalFromStruct(
    const ElectronTemperatureData& shared_data) const
{
    // Convert Te from K to eV: @f$ T_{e,\text{eV}} = T_e k_B / e ==
    // T_e \cdot \text{K\_to\_eV} @f$ (units.h mirrors units.py's
    // `K_to_eV` exactly).
    const double Te_eV = shared_data.electronTemp * rizer::units::K_to_eV;
    // @f$ k(T_e) = A / ( \sqrt{T_{e,\text{eV}}} (1 + 0.27
    // T_{e,\text{eV}}^{0.55}) ) @f$; the bracketed factor is dimensionless,
    // so the result carries the same units as A (a Cantera rate
    // coefficient in kmol/m/s units).
    return m_A / (std::sqrt(Te_eV) * (1.0 + 0.27 * std::pow(Te_eV, 0.55)));
}
 
// ---------------------------------------------------------------------------------- //
//   janev-dissociative-recombination-C2Hy - protected methods                        //
// ---------------------------------------------------------------------------------- //
 
void JanevDissociativeRecombinationC2HyRate::getParameters(AnyMap& node) const
{
    if (std::isnan(m_A)) {
        return;
    }
    node["A"].setQuantity(m_A, conversionUnits());
}
 
// ---------------------------------------------------------------------------------- //
//   Druyvesteyn - constructor                                                        //
// ---------------------------------------------------------------------------------- //
 
DruyvesteynRate::DruyvesteynRate(const Cantera::AnyMap& node,
                                  const Cantera::UnitStack& rate_units)
{
    setParameters(node, rate_units);
}
 
// ---------------------------------------------------------------------------------- //
//   Druyvesteyn - public methods                                                     //
// ---------------------------------------------------------------------------------- //
 
std::unique_ptr<MultiRateBase> DruyvesteynRate::newMultiRate() const
{
    return std::make_unique<MultiRate<DruyvesteynRate, ElectronTemperatureData>>();
}
 
void DruyvesteynRate::setParameters(const AnyMap& node, const UnitStack& rate_units)
{
    ReactionRate::setParameters(node, rate_units);  // sets units + m_input
    if (!node.hasKey("A")) {
        m_A = m_b = m_Ea_K = NAN;
        return;
    }
    // A carries the rate-constant units (Te^b's contribution is treated as
    // dimensionless, exactly like the Python reference's convert_rate_coeff);
    // Ea is expressed in Kelvin so (Ea/Te)^2 is dimensionless.
    m_A = node.units().convertRateCoeff(node["A"], conversionUnits());
    m_b = node["b"].asDouble();
    m_Ea_K = node.units().convertActivationEnergy(node["Ea"], "K");
    m_valid = true;
}
 
void DruyvesteynRate::validate(const std::string& equation,
                                const Cantera::Kinetics& /*kin*/)
{
    // Mirrors the Python rate's validate. Throw std::, not
    // CanteraError("...{}...") — the fmt instantiation rule (CLAUDE.md).
    // NaN (missing "A" in the YAML rate node, see setParameters) is caught
    // separately: NaN < 0 is false, so the negative-A check below cannot
    // catch it on its own.
    if (std::isnan(m_A)) {
        throw std::domain_error(
            "DruyvesteynRate: missing pre-exponential 'A' for reaction '"
            + equation + "'.");
    }
    if (m_A < 0) {
        throw std::domain_error(
            "DruyvesteynRate: negative pre-exponential 'A' for reaction '"
            + equation + "'.");
    }
}
 
// Compiled counterpart of the Python `DruyvesteynRate` ct.ExtensibleRate
// (rizer/kin/extensible_rate.py): a modified Arrhenius in the electron
// temperature whose exponential uses the squared temperature ratio of a
// Druyvesteyn EEDF. Depends on Te only (no Tg, no reduced field).
double DruyvesteynRate::evalFromStruct(
    const ElectronTemperatureData& shared_data) const
{
    const double Te = shared_data.electronTemp;
    const double ratio = m_Ea_K / Te;
    // @f$ k(T_e) = A \, T_e^{b} \, e^{-(E_a/T_e)^2} @f$; the exponential and
    // the power factor are dimensionless, so the result carries A's units (a
    // Cantera rate coefficient in kmol/m/s units).
    return m_A * std::pow(Te, m_b) * std::exp(-ratio * ratio);
}
 
double DruyvesteynRate::logRate(double Te) const
{
    const double ratio = m_Ea_K / Te;
    return std::log(m_A) + m_b * std::log(Te) - ratio * ratio;
}
 
// ---------------------------------------------------------------------------------- //
//   Druyvesteyn - protected methods                                                  //
// ---------------------------------------------------------------------------------- //
 
void DruyvesteynRate::getParameters(AnyMap& node) const
{
    if (std::isnan(m_A)) {
        return;
    }
    node["A"].setQuantity(m_A, conversionUnits());
    node["b"] = m_b;
    node["Ea"].setQuantity(m_Ea_K, "K", true);
}
 
// ---------------------------------------------------------------------------------- //
//   reverse-two-temperature-plasma - shared data                                     //
// ---------------------------------------------------------------------------------- //
 
bool ReverseTwoTempData::update(const ThermoPhase& phase, const Kinetics& kin)
{
    if (!shadow) {
        buildShadow(phase, kin);
    }
 
    auto plasma = dynamic_cast<const PlasmaPhase*>(&phase);
    const double Tg = phase.temperature();
    const double Te = plasma ? plasma->electronTemperature() : Tg;
    const double P = phase.pressure();
    if (Tg == m_Tg && Te == m_Te && P == m_P) {
        return false;
    }
    // The shadow evaluates every heavy species' thermo fit at Tg, and at Te when
    // a reverse reaction is Te-tagged, to form K_eq, so both are held to the
    // ceiling: above it the fits extrapolate (the CH4 mechanism's C+ fits give
    // cp/R < 1 from 1.52e5 K) and the detailed-balance rates would rest on Gibbs
    // energies they cannot give.
    if (!(Tg <= TemperatureCeiling) || (anyTeTagged && !(Te <= TemperatureCeiling))) {
        throw TrialStateError("ReverseTwoTempData::update",
            "Tg = " + std::to_string(Tg) + " K, Te = " + std::to_string(Te) +
            " K; detailed-balance reverse rates evaluate the heavy species' "
            "thermo fits at both, which hold only up to " +
            std::to_string(TemperatureCeiling) + " K.");
    }
    m_Tg = Tg;
    m_Te = Te;
    m_P = P;
 
    // Push the composition onto the ideal-gas shadow, then read the forward rate
    // constants and (concentration) equilibrium constants at BOTH Tg and Te. A
    // two-temperature-plasma forward rate in this single-T phase sees Tg=Te=T, so
    // evaluating at Te gives kf(Te) and at Tg gives kf(Tg) -- exactly what the
    // Python reference's ideal-gas shadow does (gas.TP = T, P).
    auto sth = shadow->thermo();
    auto skin = shadow->kinetics();
    std::vector<double> Y(phase.nSpecies());
    for (size_t k = 0; k < phase.nSpecies(); k++) {
        Y[k] = phase.massFraction(k);
    }
    sth->setMassFractions(Y);
    sth->setState_TP(Tg, P);
    skin->getFwdRateConstants(kf_Tg);
    skin->getEquilibriumConstants(Keq_Tg);
    // ln K is read only where kf or K is not a normal number (the log-space
    // path of evalFromStruct), so it is formed only when that can happen.
    if (anyNonNormal(kf_Tg, Keq_Tg)) {
        computeLogEquilibriumConstants(logKeq_Tg);
    }
    if (!anyTeTagged) {
        return true;
    }
    sth->setState_TP(Te, P);
    skin->getFwdRateConstants(kf_Te);
    skin->getEquilibriumConstants(Keq_Te);
    if (anyNonNormal(kf_Te, Keq_Te)) {
        computeLogEquilibriumConstants(logKeq_Te);
    }
    return true;
}
 
void ReverseTwoTempData::buildShadow(const ThermoPhase& phase, const Kinetics& kin)
{
    // Thermo: a single-temperature IDEAL-GAS copy of the main phase (electron is
    // a normal ideal-gas species), mirroring rizer's Python shadow
    // (ct.Solution(thermo="ideal-gas", species=plasma.species())). The reverse
    // rate uses detailed balance at one temperature (Tg or Te), NOT a 2-T K_eq.
    auto thermo = std::make_shared<IdealGasPhase>();
    // Auto-add elements (incl. "E") as species are added.
    thermo->addUndefinedElements();
    for (size_t k = 0; k < phase.nSpecies(); k++) {
        thermo->addSpecies(phase.species(k));
    }
    thermo->initThermo();
 
    // Kinetics: every "forward" reaction (i.e. everything EXCEPT the custom
    // reverse-two-temperature-plasma rates). We keep the janev and
    // two-temperature-plasma forwards because some reverse reactions name a
    // janev forward in their `forward_equation`, so its k_fwd/K_eq must be
    // computable here. Excluding only the reverse rates is what stops the
    // shadow from recursing back into this code.
    auto skin = newKinetics("bulk");
    skin->addThermo(thermo);
    // Size species arrays from the phase (factory does this before
    // addReaction).
    skin->init();
    // Reaction::rate() is non-const; we only read reactions, so take a non-const
    // view of the (const) source kinetics. The reaction objects are SHARED with
    // the main kinetics (as rizer's Python shadow phase does): the native rates'
    // setContext is stateless (TwoTempPlasmaRate only checks reversibility) and
    // MultiRate evaluators key off their own internal indices, not the shared
    // rate's index, so sharing does not corrupt the main solution. Sharing also
    // avoids a parameters()->newReaction round-trip, which fails because A is
    // serialized as a units-bearing Quantity that convert() cannot re-read.
    auto& mkin = const_cast<Kinetics&>(kin);
    for (size_t i = 0; i < mkin.nReactions(); i++) {
        auto r = mkin.reaction(i);
        std::string t = r->rate() ? r->rate()->type() : std::string();
        if (t == "reverse-two-temperature-plasma") {
            if (r->rate()->parameters().getString("T", "Tg") == "Te") {
                anyTeTagged = true;
            }
            continue;
        }
        skin->addReaction(r, false);
        // A second forward reaction with the same equation makes the lookup
        // ambiguous; the index is withdrawn so a reverse rate naming it raises.
        if (fwdIndex.count(r->equation()) != 0) {
            ambiguousForward.insert(r->equation());
        }
        fwdIndex[r->equation()] = skin->nReactions() - 1;
        fwdRate.push_back(r->rate());
    }
    skin->resizeReactions();
 
    shadow = Solution::create();
    shadow->setThermo(thermo);
    shadow->setKinetics(skin);
    const std::size_t nR = skin->nReactions();
    kf_Tg.assign(nR, 0.0);  Keq_Tg.assign(nR, 0.0);
    kf_Te.assign(nR, 0.0);  Keq_Te.assign(nR, 0.0);
    logKeq_Tg.assign(nR, 0.0);  logKeq_Te.assign(nR, 0.0);
    m_deltaG.assign(nR, 0.0);
    deltaN.assign(nR, 0.0);
    for (std::size_t i = 0; i < nR; i++) {
        for (std::size_t k = 0; k < thermo->nSpecies(); k++) {
            deltaN[i] += skin->productStoichCoeff(k, i) - skin->reactantStoichCoeff(k, i);
        }
    }
}
 
void ReverseTwoTempData::computeLogEquilibriumConstants(std::vector<double>& logKeq)
{
    auto sth = shadow->thermo();
    shadow->kinetics()->getDeltaSSGibbs(m_deltaG);
    const double logStandardConcentration = std::log(sth->standardConcentration());
    for (std::size_t i = 0; i < logKeq.size(); i++) {
        logKeq[i] = -m_deltaG[i] / sth->RT() + deltaN[i] * logStandardConcentration;
    }
}
 
// ---------------------------------------------------------------------------------- //
//   reverse-two-temperature-plasma - constructor                                     //
// ---------------------------------------------------------------------------------- //
 
ReverseTwoTemperaturePlasmaRate::ReverseTwoTemperaturePlasmaRate(
    const Cantera::AnyMap& node, const Cantera::UnitStack& rate_units)
{
    setParameters(node, rate_units);
}
 
// ---------------------------------------------------------------------------------- //
//   reverse-two-temperature-plasma - public methods                                  //
// ---------------------------------------------------------------------------------- //
 
std::unique_ptr<MultiRateBase>
ReverseTwoTemperaturePlasmaRate::newMultiRate() const
{
    return std::make_unique<
        MultiRate<ReverseTwoTemperaturePlasmaRate, ReverseTwoTempData>>();
}
 
void ReverseTwoTemperaturePlasmaRate::setParameters(
    const AnyMap& node, const UnitStack& rate_units)
{
    ReactionRate::setParameters(node, rate_units);
    m_forward_equation = node.getString("forward_equation", "");
    m_T = node.getString("T", "Tg");
    m_idx = npos;
    m_valid = !m_forward_equation.empty();
}
 
// Compiled counterpart of the Python `ReverseTwoTemperaturePlasma`
// ct.ExtensibleRate (rizer/kin/extensible_rate.py). Implements detailed
// balance at a single (Tg- or Te-) temperature for the reverse of a forward
// two-temperature-plasma-family reaction:
// @f[
//     k_{rev}(T) = \frac{k_{fwd}(T)}{K_{eq}(T)}
// @f]
// with T = Te for electron-impact reactions and T = Tg for heavy-species
// reactions (per the reaction's `T:` tag). k_fwd and K_eq are both looked up
// from the shadow ideal-gas solution built/cached in ReverseTwoTempData.
double ReverseTwoTemperaturePlasmaRate::evalFromStruct(
    const ReverseTwoTempData& shared_data) const
{
    if (m_idx == npos) {
        auto it = shared_data.fwdIndex.find(m_forward_equation);
        if (it == shared_data.fwdIndex.end()) {
            throw CanteraError("ReverseTwoTemperaturePlasmaRate::evalFromStruct",
                "Forward reaction '{}' not found in the forward (shadow) mechanism.",
                m_forward_equation);
        }
        if (shared_data.ambiguousForward.count(m_forward_equation) != 0) {
            throw std::invalid_argument(
                "ReverseTwoTemperaturePlasmaRate: forward_equation '" +
                m_forward_equation + "' matches more than one forward reaction "
                "(duplicates), so the reverse rate has no unique forward rate.");
        }
        m_idx = it->second;
    }
    // Detailed balance at the tagged temperature:
    // @f$ k_{rev} = k_f(T) / K_c(T) @f$.
    const bool useTe = (m_T == "Te");
    const std::vector<double>& kf = useTe ? shared_data.kf_Te : shared_data.kf_Tg;
    const std::vector<double>& Keq = useTe ? shared_data.Keq_Te : shared_data.Keq_Tg;
    const double k = kf[m_idx];
    const double K = Keq[m_idx];
    const double T = useTe ? shared_data.electronTemperature()
                           : shared_data.gasTemperature();
    double kr;
    if (std::isnormal(k) && std::isnormal(K)) {
        kr = k / K;
    } else {
        // K underflows to 0 for the electron-impact ionizations of the CH4
        // mechanism below about 746 K (C++ -> C+++), and kf underflows with it
        // (to 0, or to a subnormal within a few K of where K does), while their
        // ratio, the three-body recombination rate, is finite. Take it in log
        // space.
        const std::vector<double>& logKeq =
            useTe ? shared_data.logKeq_Te : shared_data.logKeq_Tg;
        const double logkf = logForwardRateConstant(*shared_data.fwdRate[m_idx], k,
                                                    T, m_forward_equation);
        kr = std::exp(logkf - logKeq[m_idx]);
    }
    if (!std::isfinite(kr)) {
        throw TrialStateError("ReverseTwoTemperaturePlasmaRate::evalFromStruct",
            "the reverse rate of '" + m_forward_equation + "' overflows at T = " +
            std::to_string(T) + " K.");
    }
    return kr;
}
 
// ---------------------------------------------------------------------------------- //
//   reverse-two-temperature-plasma - protected methods                               //
// ---------------------------------------------------------------------------------- //
 
void ReverseTwoTemperaturePlasmaRate::getParameters(AnyMap& node) const
{
    node["forward_equation"] = m_forward_equation;
    node["T"] = m_T;
}
 
} // namespace rizer