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 "cantera/kinetics/ReactionRateFactory.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>
 
using namespace Cantera;
 
namespace rizer {
 
// ===========================================================================
//  Registration
// ===========================================================================
 
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
// ===========================================================================
 
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;
}
 
JanevDissociativeRecombinationC2HyRate::JanevDissociativeRecombinationC2HyRate(
    const Cantera::AnyMap& node, const Cantera::UnitStack& rate_units)
{
    setParameters(node, rate_units);
}
 
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;
}
 
void JanevDissociativeRecombinationC2HyRate::getParameters(AnyMap& node) const
{
    if (std::isnan(m_A)) {
        return;
    }
    node["A"].setQuantity(m_A, conversionUnits());
}
 
// 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)));
}
 
 
// ===========================================================================
//  Druyvesteyn
// ===========================================================================
 
DruyvesteynRate::DruyvesteynRate(const Cantera::AnyMap& node,
                                  const Cantera::UnitStack& rate_units)
{
    setParameters(node, rate_units);
}
 
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 + "'.");
    }
}
 
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);
}
 
// 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);
}
 
 
// ===========================================================================
//  reverse-two-temperature-plasma
// ===========================================================================
 
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") {
            continue;
        }
        skin->addReaction(r, false);
        fwdIndex[r->equation()] = skin->nReactions() - 1;
    }
    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);
}
 
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;
    }
    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);
    sth->setState_TP(Te, P);
    skin->getFwdRateConstants(kf_Te);
    skin->getEquilibriumConstants(Keq_Te);
    return true;
}
 
ReverseTwoTemperaturePlasmaRate::ReverseTwoTemperaturePlasmaRate(
    const Cantera::AnyMap& node, const Cantera::UnitStack& rate_units)
{
    setParameters(node, rate_units);
}
 
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();
}
 
void ReverseTwoTemperaturePlasmaRate::getParameters(AnyMap& node) const
{
    node["forward_equation"] = m_forward_equation;
    node["T"] = m_T;
}
 
// 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);
        }
        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 = Keq[m_idx];
    return (K > 0.0) ? kf[m_idx] / K : 0.0;
}
 
} // namespace rizer