From 2a65df058f08984c29634d814201fed9b05b411b Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 14 May 2026 13:19:07 +0200 Subject: [PATCH 01/27] [kinetics] Add VibrationalRelaxationRate reaction rate Add the VibrationalRelaxationRate class for vibrational relaxation reactions, registered under the `vibrational-relaxation` rate type. It supports four models: `multi-state-resolved` (detailed V-V/V-T relaxation with an SSH harmonic-oscillator scaling of a base rate constant) and the mean-vibrational-energy V-T relaxation models `constant`, `castela`, and `starikovskiy`. See the class documentation for the governing rate expression and references. Also generalize TwoTempPlasmaRate with an optional gas-temperature exponent `b-gas` (default 0, which recovers the Kossyi form); a non-zero value applies the generalized k_f = A T^{b_g} T_e^b exp(-Ea_g/RT) exp(Ea_e(Te-T)/(R T Te)) expression. --- include/cantera/kinetics/TwoTempPlasmaRate.h | 20 +- .../kinetics/VibrationalRelaxationRate.h | 293 ++++++ src/kinetics/ReactionRateFactory.cpp | 6 + src/kinetics/TwoTempPlasmaRate.cpp | 17 +- src/kinetics/VibrationalRelaxationRate.cpp | 835 ++++++++++++++++++ test/kinetics/kineticsFromScratch.cpp | 1 + test/kinetics/kineticsFromYaml.cpp | 1 + 7 files changed, 1166 insertions(+), 7 deletions(-) create mode 100644 include/cantera/kinetics/VibrationalRelaxationRate.h create mode 100644 src/kinetics/VibrationalRelaxationRate.cpp diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 7f17e9de723..441da3bf33e 100644 --- a/include/cantera/kinetics/TwoTempPlasmaRate.h +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -48,12 +48,14 @@ struct TwoTempPlasmaData : public ReactionData * activation energy for electron is included. * * @f[ - * k_f = A T_e^b \exp (-E_{a,g}/RT) \exp (E_{a,e} (T_e - T)/(R T T_e)) + * k_f = A T^{b_g} T_e^b * exp(-E_{a,g}/RT) * exp(E_{a,e}(T_e - T)/(R T T_e)) * @f] * * where @f$ T_e @f$ is the electron temperature, @f$ E_{a,g} @f$ is the activation * energy for gas, and @f$ E_{a,e} @f$ is the activation energy for electron, see * Kossyi, et al. @cite kossyi1992. + * The optional gas temperature exponent b_g defaults to zero, which stricly corresponds to @cite kossyi1992. + * If b_g is non-zero, a generalisation is used. * * @ingroup arrheniusGroup */ @@ -66,10 +68,13 @@ class TwoTempPlasmaRate : public ArrheniusBase /*! * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units * depend on the reaction order and the dimensionality (surface or bulk). - * @param b Temperature exponent (non-dimensional) + * @param b Electron temperature exponent (non-dimensional) * @param Ea Activation energy in energy units [J/kmol] * @param EE Activation electron energy in energy units [J/kmol] + * @param bg Gas temperature exponent (non-dimensional). If not specified, defaults to 0. */ + TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg); + TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}); @@ -90,12 +95,12 @@ class TwoTempPlasmaRate : public ArrheniusBase */ double evalFromStruct(const TwoTempPlasmaData& shared_data) const { // m_E4_R is the electron activation (in temperature units) - return m_A * std::exp(m_b * shared_data.logTe - - m_Ea_R * shared_data.recipT + - m_E4_R * (shared_data.electronTemp - shared_data.temperature) - * shared_data.recipTe * shared_data.recipT); + return m_A * std::exp(m_bg * shared_data.logT + m_b * shared_data.logTe + - m_Ea_R * shared_data.recipT + m_E4_R * (shared_data.electronTemp - shared_data.temperature) + * shared_data.recipTe * shared_data.recipT); } + //! Evaluate derivative of reaction rate with respect to temperature //! divided by reaction rate /*! @@ -109,6 +114,9 @@ class TwoTempPlasmaRate : public ArrheniusBase double activationElectronEnergy() const { return m_E4_R * GasConstant; } + +protected: + double m_bg = 0.0; //!< Gas temperature exponent }; } diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h new file mode 100644 index 00000000000..2df4def6e8a --- /dev/null +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -0,0 +1,293 @@ +//! @file VibrationalRelaxationRate.h +//! Header for vibrational relaxation reaction rates in plasma kinetics. +//! +//! This file is part of Cantera. See License.txt in the top-level directory or +//! at https://cantera.org/license.txt for license and copyright information. + + +//! It implements the relaxation rate of vibrationally excited species in plasma kinetics. +//! Four models are currently supported by this class: 'constant', 'multi-state-resolved', 'starikovskiy', and 'castela'. + +//! The 'constant' model relaxes the vibrational species with a constant rate coefficient. It could just as well be an +//! Arrhenius rate, but the constant model is provided for convenience and to avoid confusion with conventional Arrhenius rates. + +//! The 'multi-state-resolved' model is meant to fully resolve vibrational relaxation by taking into account all vibrational species +//! in the phase (Ex: N2(v=1-8)) and solves for their V-T and V-V relaxation. The scaling of the rates are based on the SSH theory +//! detailed in Chapter 7 of the following book: +//! Capitelli, M., Ferreira, C. M., Gordiets, B. F., & Osipov, A. I. (2013). +//! Plasma kinetics in atmospheric gases (Vol. 31). Springer Science & Business Media. +//! The simplified version of the SSH theory implemented in this model is based on the harmonic oscillator approximation and can be found +//! in the following paper (equations 18 and 19): +//! Vasco Guerra et al 2019 Plasma Sources Sci. Technol. 28 073001 +//! The k10 rates are taken from //! Zhong, H. et al. (2023) “Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted +//! low-temperature NH3 oxidation,” Combustion and Flame, 256, which itself took those rates from both the Capitelli book and +//! the article from Strikovskiy and Aleksandrov (2013), see more below. + +//! The 'castela' model is meant to be used only for N2 vibrational relaxation, by collisions with N2, O2 and O exclusively. +//! It implements the mean vibrational energy relaxation model using a fictive cantera species. +//! The Castela model is based on the following paper: +//! Castela, M., Fiorina, B., Coussement, A., Gicquel, O., Darabiha, N., & Laux, C. O. (2016). +//! Modelling the impact of non-equilibrium discharges on reactive mixtures for simulations of +//! plasma-assisted ignition in turbulent flows. Combustion and flame, 166, 133-147. + +//! The so-called 'starikovskiy' model is an extension of the Castela model to several vibrational species and additional colliders. +//! A lot of vibrational relaxation rates can be found in Table 1 of the following paper: +//! Starikovskiy, A., & Aleksandrov, N. (2013). Plasma-assisted ignition and combustion. +//! Progress in Energy and Combustion Science, 39(1), 61-110. +//! The rates for the vibrational relaxation of NH3 can be found in the reaction mechanism provided in supplementary material +//! of the following paper: +//! Zhong, H. et al. (2023) “Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted low-temperature NH3 oxidation,” +//! Combustion and Flame, 256. +//! More rates for the vibrational relaxation of CH$ can be found in the following paper: +//! Popov, N.A. (2016) “Kinetics of plasma-assisted combustion: Effect of non-equilibrium excitation on the ignition +//! and oxidation of combustible mixtures,” Plasma Sources Science and Technology. Institute of Physics Publishing. + + + +#ifndef CT_VIBRATIONALRELAXATIONRATE_H +#define CT_VIBRATIONALRELAXATIONRATE_H + +#include "Arrhenius.h" + +#include + +namespace Cantera +{ + +//! Shared temperature data for vibrational relaxation rates. +struct DetailedVibData : public ReactionData +{ + //! Update cached temperature-dependent data. + /*! + * @param phase Thermodynamic phase used to retrieve the gas temperature. + * @param kin Kinetics object. Not used here, but required by the + * ReactionData interface. + * + * @return `true` if the temperature has changed and rates need to be + * recomputed; `false` otherwise. + */ + bool update(const ThermoPhase& phase, const Kinetics& kin) override; + + using ReactionData::update; +}; + + +//! Vibrational relaxation reaction rate. +/*! + * This class provides a common implementation for several vibrational + * relaxation models: + * + * - `constant` + * - `multi-state-resolved` + * - `starikovskiy` + * - `castela` + * + * Internally, all models are mapped to the following generic expression: + * + * @f[ + * k_f = + * scaling \, A \, + * \exp \left( + * b \ln T + * + B + * + C T^{-1/3} + * + D T^{-m} + * + E T^{-z} + * \right) + * @f] + * + * where `T` is the gas temperature in K. + * + * The YAML reaction type is: + * + * @code{.yaml} + * type: vibrational-relaxation + * @endcode + * + * The selected physical model is specified separately using: + * + * @code{.yaml} + * vibration_model: constant + * vibration_model: multi-state-resolved + * vibration_model: starikovskiy + * vibration_model: castela + * @endcode + * + * Unit conventions: + * + * - `A` uses standard Cantera rate coefficient units. Its units depend on the + * reaction order and are converted by `ArrheniusBase`. + * - `b`, `B`, `m`, `z`, and `scaling` are dimensionless. + * - `C` is interpreted as K^(1/3), assuming `T` is in K. + * - `D` is interpreted as K^m, assuming `T` is in K. + * - `E` is interpreted as K^z, assuming `T` is in K. + * + * The coefficients `B`, `C`, `D`, `E`, `m`, `z`, and `scaling` are read as + * raw floating-point values. They are not converted by Cantera's unit system. + * + * @ingroup arrheniusGroup + */ +class VibrationalRelaxationRate : public ArrheniusBase +{ +public: + //! Default constructor. + VibrationalRelaxationRate(); + + //! Constructor using the internal generic representation. + /*! + * @param A Pre-exponential factor. + * @param B Dimensionless constant in the exponential. + * @param C Coefficient multiplying T^(-1/3). + * @param D Coefficient multiplying T^(-m). + * @param b Dimensionless temperature exponent. + * @param scaling Dimensionless scaling factor. + * @param m Temperature exponent used by the D term. + * @param E Coefficient multiplying T^(-z). + * @param z Temperature exponent used by the E term. + */ + VibrationalRelaxationRate(double A, double B, double C, double D, + double b, double scaling = 1.0, + double m = 2.0 / 3.0, + double E = 0.0, double z = 1.0); + + //! Constructor based on AnyMap content. + explicit VibrationalRelaxationRate(const AnyMap& node, + const UnitStack& rate_units = {}); + + //! Set rate parameters from an AnyMap. + void setParameters(const AnyMap& node, + const UnitStack& rate_units) override; + + //! Get rate parameters for YAML serialization. + void getParameters(AnyMap& node) const override; + + //! Create a rate evaluator for this reaction rate type. + unique_ptr newMultiRate() const override { + return make_unique>(); + } + + //! String identifying this reaction rate type. + const string type() const override { + return "vibrational-relaxation"; + } + + //! Set context of reaction rate evaluation. + /*! + * Vibrational relaxation rates are intended for irreversible + * non-equilibrium plasma reactions. Reversible reactions are rejected + * because the reverse rate cannot be obtained from conventional + * thermochemistry for these models. + */ + void setContext(const Reaction& rxn, const Kinetics& kin) override; + + //! Evaluate the forward rate coefficient. + double evalFromStruct(const DetailedVibData& shared_data) const { + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + return m_scaling * m_A * std::exp( + m_b * shared_data.logT + + m_B + + m_C * invT13 + + m_D * std::pow(invT, m_m) + + m_E * std::pow(invT, m_z) + ); + } + + //! Evaluate the scaled temperature derivative. + /*! + * This returns: + * + * @f[ + * \frac{1}{k_f} \frac{d k_f}{dT} + * = + * \frac{d \ln k_f}{dT} + * @f] + * + * For the internal generic expression, this is: + * + * @f[ + * \frac{b}{T} + * - \frac{C}{3} T^{-4/3} + * - m D T^{-m-1} + * - z E T^{-z-1} + * @f] + */ + double ddTScaledFromStruct(const DetailedVibData& shared_data) const; + +private: + //! Dimensionless constant in the exponential. + double m_B = 0.0; + + //! Coefficient multiplying T^(-1/3). + double m_C = 0.0; + + //! Coefficient multiplying T^(-m). + double m_D = 0.0; + + //! Dimensionless scaling factor. + double m_scaling = 1.0; + + //! Temperature exponent used by the D term. + double m_m = 2.0 / 3.0; + + //! Coefficient multiplying T^(-z). + double m_E = 0.0; + + //! Temperature exponent used by the E term. + double m_z = 1.0; + + //! Castela coefficient a. + double m_castela_a = 0.0; + + //! Castela coefficient b. + double m_castela_b = 0.0; + + //! Castela reference pressure. + double m_referencePressure = OneAtm; + + //! Selected vibrational relaxation model. + /*! + * Accepted values: + * + * - `constant` + * - `multi-state-resolved` + * - `starikovskiy` + * - `castela` + */ + string m_vibration_model = "multi-state-resolved"; + + //! YAML key names. + string m_B_str = "B"; + string m_C_str = "C"; + string m_D_str = "D"; + string m_scaling_str = "scaling"; + string m_m_str = "m"; + string m_E_str = "E"; + string m_z_str = "z"; + string m_reference_pressure_str = "reference-pressure"; + + //! Configure the ArrheniusBase part from an already-converted internal A value. + /*! + * This is needed for models such as Castela, where the user-facing YAML does + * not contain a standard Arrhenius A coefficient. + */ + void configureBaseFromInternalA(const AnyMap& node, + const UnitStack& rate_units, + double A, double b); + + //! Configure the ArrheniusBase part from a YAML A value and an explicit b. + /*! + * This is needed for models such as constant and Starikovskiy, where the YAML + * does not contain the standard Arrhenius pair A / b. + */ + void configureBaseFromYamlA(const AnyMap& node, + const UnitStack& rate_units, + const AnyValue& A, + double b); + + }; + +} + +#endif \ No newline at end of file diff --git a/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 1019fe4f4ed..cb079dc0ee1 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -18,6 +18,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" namespace Cantera { @@ -40,6 +41,11 @@ ReactionRateFactory::ReactionRateFactory() return new TwoTempPlasmaRate(node, rate_units); }); + // VibrationalRelaxationRate evaluator + reg("vibrational-relaxation", [](const AnyMap& node, const UnitStack& rate_units) { + return new VibrationalRelaxationRate(node, rate_units); + }); + // ElectronCollisionPlasmaRate evaluator reg("electron-collision-plasma", [](const AnyMap& node, const UnitStack& rate_units) { return new ElectronCollisionPlasmaRate(node, rate_units); diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index ce3f4694969..0b6563e98e3 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -57,19 +57,34 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) m_Ea_str = "Ea-gas"; m_E4_str = "Ea-electron"; m_E4_R = EE / GasConstant; + m_bg = 0.0; +} + +TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg) + : ArrheniusBase(A, b, Ea) +{ + TwoTempPlasmaRate(A, b, Ea, EE); + m_bg = bg; } TwoTempPlasmaRate::TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units) : TwoTempPlasmaRate() { setParameters(node, rate_units); + + if (node.hasKey("b-gas")) { + m_bg = node["b-gas"].asDouble(); + } else if (node.hasKey("b_gas")) { + m_bg = node["b_gas"].asDouble(); + } } double TwoTempPlasmaRate::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const { warn_user("TwoTempPlasmaRate::ddTScaledFromStruct", "Temperature derivative does not consider changes of electron temperature."); - return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; + return m_bg * shared_data.recipT + + (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; } void TwoTempPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp new file mode 100644 index 00000000000..79aabf911ee --- /dev/null +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -0,0 +1,835 @@ +//! @file VibrationalRelaxationRate.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/kinetics/VibrationalRelaxationRate.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/base/global.h" + +#include +#include +#include +#include +#include +#include + +namespace Cantera +{ + +namespace +{ + +void requireNoKey(const AnyMap& node, const string& key, + const string& model, const string& where) +{ + if (node.hasKey(key)) { + throw InputFileError(where, node, + "Key '{}' is not allowed for vibration_model '{}'.", key, model); + } +} + + +void requireKey(const AnyMap& node, const string& key, + const string& model, const string& where) +{ + if (!node.hasKey(key)) { + throw InputFileError(where, node, + "Missing required key '{}' for vibration_model '{}'.", key, model); + } +} + + +bool isVibrationalSpecies(const string& name) +{ + // Supported names: + // + // N2(v) + // O2(v0) + // O2(v1) + // O2(v12) + // + // The old form O2(v=1) is intentionally not required. + const auto pos = name.find("(v"); + return pos != string::npos && !name.empty() && name.back() == ')'; +} + + +string groundStateName(const string& name) +{ + const auto pos = name.find("(v"); + if (pos == string::npos) { + return name; + } + return name.substr(0, pos); +} + + +string vibrationalFamilyName(const string& name) +{ + // Collapse O2(v1), O2(v2), O2(v12), and O2(v) into O2(v). + return groundStateName(name) + "(v)"; +} + + +double compositionSum(const Composition& comp) +{ + double sum = 0.0; + for (const auto& item : comp) { + sum += item.second; + } + return sum; +} + + +Composition replaceVibrationalSpeciesByGroundState(const Composition& comp) +{ + Composition out; + for (const auto& item : comp) { + const string& name = item.first; + const double value = item.second; + + if (isVibrationalSpecies(name)) { + out[groundStateName(name)] += value; + } else { + out[name] += value; + } + } + return out; +} + + +bool sameComposition(const Composition& a, const Composition& b, + double tol = 1e-12) +{ + std::set names; + + for (const auto& item : a) { + names.insert(item.first); + } + for (const auto& item : b) { + names.insert(item.first); + } + + for (const auto& name : names) { + double av = 0.0; + double bv = 0.0; + + const auto ait = a.find(name); + if (ait != a.end()) { + av = ait->second; + } + + const auto bit = b.find(name); + if (bit != b.end()) { + bv = bit->second; + } + + if (std::abs(av - bv) > tol) { + return false; + } + } + + return true; +} + + +std::vector vibrationalSpeciesInComposition(const Composition& comp) +{ + std::vector out; + + for (const auto& item : comp) { + const string& name = item.first; + const double value = item.second; + + if (value != 0.0 && isVibrationalSpecies(name)) { + out.push_back(name); + } + } + + return out; +} + + +// Registry used to check that a given vibrational family uses exactly one +// relaxation model inside one Kinetics object. +// +// Allowed: +// +// N2(v) -> constant +// O2(v) -> multi-state-resolved +// NH3(v) -> starikovskiy +// +// Forbidden: +// +// N2(v) + O -> castela +// N2(v) + N2 -> starikovskiy +// +std::map> s_modelByFamily; +std::set s_warnedMixedModels; + + +void registerVibrationalModelConsistency(const Kinetics& kin, + const string& family, + const string& model, + const AnyMap& input) +{ + auto& modelByFamily = s_modelByFamily[&kin]; + + const auto existing = modelByFamily.find(family); + if (existing != modelByFamily.end() && existing->second != model) { + throw InputFileError("VibrationalRelaxationRate::setContext", input, + "Inconsistent vibration_model for vibrational family '{}'. " + "This family was already registered with model '{}', but the " + "current reaction uses model '{}'. A given vibrational family " + "must use exactly one relaxation model.", + family, existing->second, model); + } + + modelByFamily[family] = model; + + std::set models; + for (const auto& item : modelByFamily) { + models.insert(item.second); + } + + if (models.size() > 1 && !s_warnedMixedModels.count(&kin)) { + std::ostringstream msg; + msg << "Multiple vibrational relaxation models were detected in the " + << "same kinetics object. This is allowed only if each " + << "vibrational family is internally consistent:\n"; + + for (const auto& item : modelByFamily) { + msg << " - " << item.first << ": " << item.second << "\n"; + } + + warn_user("VibrationalRelaxationRate::setContext", msg.str()); + s_warnedMixedModels.insert(&kin); + } +} + + +string inferRelaxingFamily(const Reaction& rxn) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + + if (vibReactants.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "A vibrational-relaxation reaction must contain at least one " + "vibrational reactant, for example O2(v1), N2(v), or NH3(v1)."); + } + + return vibrationalFamilyName(vibReactants.front()); +} + + +void validateSimpleRelaxationToGroundState(const Reaction& rxn, + const string& model) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); + + if (vibReactants.size() != 1) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' expects exactly one vibrational reactant.", + model); + } + + if (!vibProducts.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' describes relaxation to the ground state " + "and therefore does not allow vibrational products.", model); + } + + const Composition relaxedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + if (!sameComposition(relaxedReactants, rxn.products)) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid vibrational relaxation stoichiometry for model '{}'. " + "Expected a reaction equivalent to X(v) + M => X + M, where the " + "collider M is unchanged.", model); + } + + if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 + || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) + { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' expects a bimolecular relaxation reaction " + "of the form X(v) + M => X + M.", model); + } +} + + +void validateCastelaReaction(const Reaction& rxn) +{ + validateSimpleRelaxationToGroundState(rxn, "castela"); + + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const string family = vibrationalFamilyName(vibReactants.front()); + + if (family != "N2(v)") { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "The Castela relaxation model is only valid for N2 vibrational " + "relaxation. Found vibrational family '{}'.", family); + } + + // Castela parameters used here are only intended for the following + // colliders: N2, O2, and O. + const Composition collapsedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + string collider = ""; + for (const auto& item : collapsedReactants) { + const string& name = item.first; + const double value = item.second; + + if (name == "N2") { + if (std::abs(value - 2.0) < 1e-12) { + collider = "N2"; + } + } else if (std::abs(value - 1.0) < 1e-12) { + collider = name; + } + } + + if (collider != "N2" && collider != "O2" && collider != "O") { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "The Castela relaxation model only supports colliders 'N2', " + "'O2', and 'O'. Found collider '{}'.", collider); + } +} + + +void validateDetailedRelaxationReaction(const Reaction& rxn) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); + + if (vibReactants.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "A detailed vibrational relaxation reaction must contain at least " + "one vibrational reactant."); + } + + // All vibrational species in a detailed relaxation reaction are required + // to belong to the same vibrational family. + const string family = vibrationalFamilyName(vibReactants.front()); + + for (const auto& sp : vibReactants) { + const string spFamily = vibrationalFamilyName(sp); + if (spFamily != family) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: all " + "vibrational reactants must belong to the same vibrational " + "family. Found '{}' and '{}'.", family, spFamily); + } + } + + for (const auto& sp : vibProducts) { + const string spFamily = vibrationalFamilyName(sp); + if (spFamily != family) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: all " + "vibrational products must belong to the same vibrational " + "family. Found '{}' and '{}'.", family, spFamily); + } + } + + // The reaction must conserve the ground-state composition once all + // vibrational labels are removed. + const Composition collapsedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + const Composition collapsedProducts = + replaceVibrationalSpeciesByGroundState(rxn.products); + + if (!sameComposition(collapsedReactants, collapsedProducts)) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: replacing all " + "vibrational species by their ground-state species does not " + "conserve stoichiometry."); + } + + // The current detailed VV/VT formulation is bimolecular. + if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 + || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) + { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: expected a " + "bimolecular reaction with two reactant molecules and two product " + "molecules."); + } +} + +} // namespace + + +bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + const double T = phase.temperature(); + + if (T == temperature) { + return false; + } + + ReactionData::update(T); + return true; +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate() +{ +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate( + double A, double B, double C, double D, + double b, double scaling, double m, double E, double z) + : ArrheniusBase(A, b, 0.0) + , m_B(B) + , m_C(C) + , m_D(D) + , m_scaling(scaling) + , m_m(m) + , m_E(E) + , m_z(z) +{ +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate( + const AnyMap& node, const UnitStack& rate_units) + : VibrationalRelaxationRate() +{ + setParameters(node, rate_units); +} + + +void VibrationalRelaxationRate::configureBaseFromInternalA( + const AnyMap& node, const UnitStack& rate_units, double A, double b) +{ + // Store the original input and configure reaction-rate units. + // + // We intentionally do not call ArrheniusBase::setParameters here because + // some vibration models do not expose a standard YAML rate-constant with + // both A and b. Castela is the main example. + ReactionRate::setParameters(node, rate_units); + setRateUnits(rate_units); + + m_negativeA_ok = node.getBool("negative-A", false); + + m_A = A; + m_b = b; + + if (m_A != 0.0) { + m_logA = std::log(std::abs(m_A)); + } else { + m_logA = NAN; + } + + // Critical: ArrheniusBase::validate checks this flag. + m_valid = true; +} + + +void VibrationalRelaxationRate::configureBaseFromYamlA( + const AnyMap& node, const UnitStack& rate_units, + const AnyValue& A, double b) +{ + // Store the original input and configure reaction-rate units first, so + // conversionUnits() is available for A. + ReactionRate::setParameters(node, rate_units); + setRateUnits(rate_units); + + m_negativeA_ok = node.getBool("negative-A", false); + + // Convert the user-facing YAML A value with Cantera's standard + // rate-coefficient unit conversion. + m_A = node.units().convertRateCoeff(A, conversionUnits()); + m_b = b; + + if (m_A != 0.0) { + m_logA = std::log(std::abs(m_A)); + } else { + m_logA = NAN; + } + + // Critical: ArrheniusBase::validate checks this flag. + m_valid = true; +} + + +void VibrationalRelaxationRate::setParameters(const AnyMap& node, + const UnitStack& rate_units) +{ + // The reaction rate type is always: + // + // type: vibrational-relaxation + // + // The physical model is selected separately by: + // + // vibration_model: constant + // vibration_model: multi-state-resolved + // vibration_model: starikovskiy + // vibration_model: castela + // + // For backward compatibility, the default model is multi-state-resolved. + m_vibration_model = "multi-state-resolved"; + + if (node.hasKey("vibration_model")) { + m_vibration_model = node["vibration_model"].asString(); + } else if (node.hasKey("model")) { + m_vibration_model = node["model"].asString(); + } + + if (!node.hasKey("rate-constant")) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "A vibrational-relaxation reaction requires a 'rate-constant' " + "mapping."); + } + + const auto& rate = node["rate-constant"]; + + if (!rate.is()) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "The 'rate-constant' field must be a mapping."); + } + + const auto& rateMap = rate.as(); + + if (m_vibration_model == "constant") { + // Constant model: + // + // k(T) = A + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_b_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_B_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_C_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_D_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); + + m_B = 0.0; + m_C = 0.0; + m_D = 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = 1.0; + } + else if (m_vibration_model == "multi-state-resolved") { + // Detailed VV/VT model: + // + // k(T) = scaling * A * exp( + // b * log(T) + // + B + // + C * T^(-1/3) + // + D * T^(-2/3) + // ) + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + ArrheniusBase::setParameters(node, rate_units); + + m_B = rateMap.hasKey(m_B_str) ? rateMap[m_B_str].asDouble() : 0.0; + m_C = rateMap.hasKey(m_C_str) ? rateMap[m_C_str].asDouble() : 0.0; + m_D = rateMap.hasKey(m_D_str) ? rateMap[m_D_str].asDouble() : 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = rateMap.hasKey(m_scaling_str) + ? rateMap[m_scaling_str].asDouble() + : 1.0; + } + else if (m_vibration_model == "starikovskiy") { + // User-facing formula: + // + // k(T) = A * T^n * exp( + // K + // - B * T^(-1/3) + // + C * T^(-m) + // + D * T^(-z) + // ) + // + // Internal mapping: + // + // m_b = n + // m_B = K + // m_C = -B + // m_D = C + // m_m = m + // m_E = D + // m_z = z + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_b_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + const double n = rateMap.hasKey("n") ? rateMap["n"].asDouble() : 0.0; + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], n); + + m_B = rateMap.hasKey("K") ? rateMap["K"].asDouble() : 0.0; + m_C = rateMap.hasKey("B") ? -rateMap["B"].asDouble() : 0.0; + m_D = rateMap.hasKey("C") ? rateMap["C"].asDouble() : 0.0; + m_m = rateMap.hasKey("m") ? rateMap["m"].asDouble() : 1.0; + m_E = rateMap.hasKey("D") ? rateMap["D"].asDouble() : 0.0; + m_z = rateMap.hasKey("z") ? rateMap["z"].asDouble() : 1.0; + m_scaling = 1.0; + + if (m_m <= 0.0 || m_z <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "The Starikovskiy exponents 'm' and 'z' must be positive."); + } + } + else if (m_vibration_model == "castela") { + // Castela model: + // + // Original relaxation time: + // + // tau_k = p0 / p_k + // * exp[a_k * (T^(-1/3) - b_k) - 18.42] + // + // Equivalent bimolecular rate coefficient: + // + // k_k(T) = R T / p0 + // * exp[18.42 + a_k b_k - a_k T^(-1/3)] + // + // Internal mapping: + // + // A = R / p0 + // b = 1 + // B = 18.42 + a_k b_k + // C = -a_k + // D = 0 + // E = 0 + // scaling = 1 + + requireKey(rateMap, "a", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireKey(rateMap, "b", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "K", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_B_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_C_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_D_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + m_castela_a = rateMap["a"].asDouble(); + m_castela_b = rateMap["b"].asDouble(); + + if (rateMap.hasKey(m_reference_pressure_str)) { + m_referencePressure = rateMap.convert(m_reference_pressure_str, "Pa"); + } else { + m_referencePressure = OneAtm; + } + + if (m_referencePressure <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "Castela reference-pressure must be positive."); + } + + configureBaseFromInternalA( + node, rate_units, GasConstant / m_referencePressure, 1.0); + + m_B = 18.42 + m_castela_a * m_castela_b; + m_C = -m_castela_a; + m_D = 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = 1.0; + } + else { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "Unrecognized vibration_model '{}'. Expected 'multi-state-resolved', " + "'starikovskiy', 'castela', or 'constant'.", + m_vibration_model); + } +} + + +void VibrationalRelaxationRate::getParameters(AnyMap& node) const +{ + if (!valid()) { + return; + } + + if (allowNegativePreExponentialFactor()) { + node["negative-A"] = true; + } + + node["vibration_model"] = m_vibration_model; + + AnyMap rateNode; + + auto storePreExponentialFactor = [&](AnyMap& target, double A) { + if (conversionUnits().factor() != 0.0) { + target[m_A_str].setQuantity(A, conversionUnits()); + } else { + target[m_A_str] = A; + target["__unconvertible__"] = true; + } + }; + + if (m_vibration_model == "constant") { + const double tol = 1e-12; + + if (std::abs(m_b) > tol + || std::abs(m_B) > tol + || std::abs(m_C) > tol + || std::abs(m_D) > tol + || std::abs(m_E) > tol) + { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'constant': the internal " + "parameters contain temperature-dependent terms."); + } + + storePreExponentialFactor(rateNode, m_scaling * m_A); + } + else if (m_vibration_model == "multi-state-resolved") { + storePreExponentialFactor(rateNode, m_A); + + rateNode[m_b_str] = m_b; + rateNode[m_B_str] = m_B; + rateNode[m_C_str] = m_C; + rateNode[m_D_str] = m_D; + rateNode[m_scaling_str] = m_scaling; + } + else if (m_vibration_model == "starikovskiy") { + storePreExponentialFactor(rateNode, m_A); + + rateNode["n"] = m_b; + rateNode["K"] = m_B; + rateNode["B"] = -m_C; + rateNode["C"] = m_D; + rateNode["m"] = m_m; + rateNode["D"] = m_E; + rateNode["z"] = m_z; + } + else if (m_vibration_model == "castela") { + const double tol = 1e-12; + + if (std::abs(m_b - 1.0) > tol + || std::abs(m_D) > tol + || std::abs(m_E) > tol + || std::abs(m_scaling - 1.0) > tol) + { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'castela': the internal " + "parameters are not consistent with the Castela form. " + "Expected b = 1, D = 0, E = 0, and scaling = 1."); + } + + if (m_referencePressure <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'castela': " + "reference-pressure must be positive."); + } + + rateNode["a"] = m_castela_a; + rateNode["b"] = m_castela_b; + rateNode[m_reference_pressure_str].setQuantity(m_referencePressure, "Pa"); + } + else { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Unrecognized vibration_model '{}'. Expected 'multi-state-resolved', " + "'starikovskiy', 'castela', or 'constant'.", + m_vibration_model); + } + + rateNode.setFlowStyle(); + node["rate-constant"] = std::move(rateNode); +} + + +double VibrationalRelaxationRate::ddTScaledFromStruct( + const DetailedVibData& shared_data) const +{ + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + return m_b * invT + - (m_C / 3.0) * invT13 * invT + - m_m * m_D * std::pow(invT, m_m) * invT + - m_z * m_E * std::pow(invT, m_z) * invT; +} + + +void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& kin) +{ + if (rxn.reversible) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Vibrational relaxation rates do not support reversible " + "reactions."); + } + + const string family = inferRelaxingFamily(rxn); + + if (m_vibration_model == "constant") { + validateSimpleRelaxationToGroundState(rxn, "constant"); + } else if (m_vibration_model == "castela") { + validateCastelaReaction(rxn); + } else if (m_vibration_model == "starikovskiy") { + validateSimpleRelaxationToGroundState(rxn, "starikovskiy"); + } else if (m_vibration_model == "multi-state-resolved") { + validateDetailedRelaxationReaction(rxn); + } else { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Unrecognized vibration_model '{}'.", m_vibration_model); + } + + registerVibrationalModelConsistency( + kin, family, m_vibration_model, rxn.input); +} + +} // namespace Cantera \ No newline at end of file diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index b4d120099db..332ab9cb111 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -14,6 +14,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" #include "cantera/base/Array.h" #include "cantera/base/stringUtils.h" diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index fbe699b3889..d03da101282 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -13,6 +13,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/base/Array.h" From 54b42384c74349c00a3e237700681d847687ad17 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 22 May 2026 23:05:05 +0200 Subject: [PATCH 02/27] rebase and resolve conflicts - [thermo] Warn about misdeclared vibrational reservoir species Detect fictive vibrational-reservoir species (named "X(v)") and warn when the corresponding base species "X" is missing, or when the reservoir mole fraction grows into a non-negligible fraction of the base species -- a sign that the chosen reservoir energy is too low and is perturbing the chemistry. --- include/cantera/thermo/PlasmaPhase.h | 15 +++++ src/thermo/PlasmaPhase.cpp | 98 ++++++++++++++++++++++++++++ 2 files changed, 113 insertions(+) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 00907cb0a90..0b3b1e2bce4 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -1025,6 +1025,21 @@ class PlasmaPhase: public IdealGasPhase //! Work array mutable std::vector m_work; + + struct VibrationalReservoirSpecies { + size_t reservoirIndex = npos; + size_t baseSpeciesIndex = npos; + }; + + std::vector m_vibrationalReservoirSpecies; + + bool m_vibrationalReservoirSpeciesNeedUpdate = true; + + double m_vibrationalMoleFractionThreshold = 1e-2; + double m_vibrationalAbsoluteMoleFractionThreshold = 1e-20; + + void updateVibrationalReservoirSpecies(); + void checkVibrationalReservoirMoleFractions(); }; } diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 48eeb01bfe0..69c6cab2cfd 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -19,6 +19,22 @@ namespace Cantera { namespace { const double gamma = sqrt(2 * ElectronCharge / ElectronMass); + + bool isVibrationalReservoirName(const string& name, string& baseName) + { + const string suffix = "(v)"; + + if (name.size() <= suffix.size()) { + return false; + } + + if (name.compare(name.size() - suffix.size(), suffix.size(), suffix) != 0) { + return false; + } + + baseName = name.substr(0, name.size() - suffix.size()); + return !baseName.empty(); + } } PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) @@ -109,6 +125,11 @@ bool PlasmaPhase::addSpecies(shared_ptr spec) "Only one electron species is allowed.", spec->name); } } + + if (added) { + m_vibrationalReservoirSpeciesNeedUpdate = true; + } + return added; } @@ -1006,6 +1027,83 @@ double PlasmaPhase::intrinsicHeating() const double qJ = jouleHeatingPower(); checkFinite(qJ); + checkVibrationalReservoirMoleFractions(); + return qJ; } + +void PlasmaPhase::updateVibrationalReservoirSpecies() +{ + m_vibrationalReservoirSpecies.clear(); + + for (size_t k = 0; k < nSpecies(); k++) { + string baseName; + const string& reservoirName = speciesName(k); + + if (!isVibrationalReservoirName(reservoirName, baseName)) { + continue; + } + + size_t kBase = speciesIndex(baseName, false); + + if (kBase == npos) { + warn_user("PlasmaPhase::updateVibrationalReservoirSpecies", + "Species '{}' matches the fictive vibrational-reservoir naming " + "convention, but the corresponding base species '{}' was not found. " + "This species will not be monitored as a vibrational reservoir.", + reservoirName, baseName); + continue; + } + + VibrationalReservoirSpecies reservoir; + reservoir.reservoirIndex = k; + reservoir.baseSpeciesIndex = kBase; + + m_vibrationalReservoirSpecies.push_back(reservoir); + } + + m_vibrationalReservoirSpeciesNeedUpdate = false; +} + +void PlasmaPhase::checkVibrationalReservoirMoleFractions() +{ + if (m_vibrationalReservoirSpeciesNeedUpdate) { + updateVibrationalReservoirSpecies(); + } + + for (const auto& reservoir : m_vibrationalReservoirSpecies) { + const size_t kReservoir = reservoir.reservoirIndex; + const size_t kBase = reservoir.baseSpeciesIndex; + + const double Xv = moleFraction(kReservoir); + const double Xb = moleFraction(kBase); + + const double pool = Xv + Xb; + + if (pool <= m_vibrationalAbsoluteMoleFractionThreshold) { + continue; + } + + const double reservoirFraction = Xv / pool; + + if (reservoirFraction > m_vibrationalMoleFractionThreshold) { + const string& reservoirName = speciesName(kReservoir); + const string& baseName = speciesName(kBase); + + writelog("Warning: fictive vibrational reservoir species '{}' contains " + "{:.3e} of the total '{}' pool. " + "X({}) = {:.3e}, X({}) = {:.3e}, threshold = {:.3e}. " + "Chemistry involving '{}' may be affected because part of the " + "material is stored in an inert vibrational reservoir.\n", + reservoirName, + reservoirFraction, + baseName, + reservoirName, Xv, + baseName, Xb, + m_vibrationalMoleFractionThreshold, + baseName); + } + } +} + } From abd1c0992308724f25749a7095fb126ed53d719e Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 12:33:59 +0200 Subject: [PATCH 03/27] Updated TwoTempPlasmaRate class comments and formatting --- include/cantera/kinetics/TwoTempPlasmaRate.h | 41 +++++++++++++------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 441da3bf33e..581e2899841 100644 --- a/include/cantera/kinetics/TwoTempPlasmaRate.h +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -12,7 +12,7 @@ namespace Cantera { -//! Data container holding shared data specific to TwoTempPlasmaRate +//! Data container holding shared data specific to TwoTempPlasmaRate. /** * The data container `TwoTempPlasmaData` holds precalculated data common to * all `TwoTempPlasmaRate` objects. @@ -38,7 +38,6 @@ struct TwoTempPlasmaData : public ReactionData double recipTe = 1.0; //!< inverse of electron temperature }; - //! Two temperature plasma reaction rate type depends on both //! gas temperature and electron temperature. /*! @@ -66,32 +65,45 @@ class TwoTempPlasmaRate : public ArrheniusBase //! Constructor. /*! - * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units + * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units. * depend on the reaction order and the dimensionality (surface or bulk). - * @param b Electron temperature exponent (non-dimensional) - * @param Ea Activation energy in energy units [J/kmol] - * @param EE Activation electron energy in energy units [J/kmol] - * @param bg Gas temperature exponent (non-dimensional). If not specified, defaults to 0. + * @param b Electron temperature exponent (non-dimensional). + * @param Ea Activation energy in energy units [J/kmol]. + * @param EE Activation electron energy in energy units [J/kmol]. + * @param bg Gas temperature exponent (non-dimensional). If not specified, defaults to 0. + * @since New in %Cantera 4.0 */ TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg); + //! Constructor. + /*! + * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units. + * depend on the reaction order and the dimensionality (surface or bulk). + * @param b Electron temperature exponent (non-dimensional). + * @param Ea Activation energy in energy units [J/kmol]. Defaults to 0. + * @param EE Activation electron energy in energy units [J/kmol]. Defaults to 0. + */ TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); + //! Constructor based on an AnyMap object instead of all parameters directly. TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}); + //! Creates a new two-temperature-plasma reaction. unique_ptr newMultiRate() const override { return make_unique>(); } + //! Returns the reaction type. const string type() const override { return "two-temperature-plasma"; } + //! Check that the reaction does not have the 'reversible: true' attribute. void setContext(const Reaction& rxn, const Kinetics& kin) override; - //! Evaluate reaction rate + //! Evaluates reaction rate. /*! - * @param shared_data data shared by all reactions of a given type + * @param shared_data data shared by all reactions of a given type. */ double evalFromStruct(const TwoTempPlasmaData& shared_data) const { // m_E4_R is the electron activation (in temperature units) @@ -100,23 +112,24 @@ class TwoTempPlasmaRate : public ArrheniusBase * shared_data.recipTe * shared_data.recipT); } - //! Evaluate derivative of reaction rate with respect to temperature - //! divided by reaction rate + //! divided by reaction rate. /*! * This method does not consider changes of electron temperature. * A corresponding warning is raised. - * @param shared_data data shared by all reactions of a given type + * @param shared_data data shared by all reactions of a given type. */ double ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const; - //! Return the electron activation energy *Ea* [J/kmol] + //! Return the electron activation energy *Ea* [J/kmol]. double activationElectronEnergy() const { return m_E4_R * GasConstant; } protected: - double m_bg = 0.0; //!< Gas temperature exponent + //! Gas temperature exponent. + //! @since New in %Cantera 4.0 + double m_bg = 0.0; }; } From 927ab86bbba63c427a1727f67016d4a07cf6a377 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 14:02:53 +0200 Subject: [PATCH 04/27] accept only b-gas as entry for yaml --- src/kinetics/TwoTempPlasmaRate.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index 0b6563e98e3..5f257fbb0a9 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -74,8 +74,6 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_u if (node.hasKey("b-gas")) { m_bg = node["b-gas"].asDouble(); - } else if (node.hasKey("b_gas")) { - m_bg = node["b_gas"].asDouble(); } } From 0ba7ea0bc1d0aacd7b772c451eae9281f227c59f Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 14:05:53 +0200 Subject: [PATCH 05/27] update two-temperature-plasma yaml documentation --- doc/sphinx/yaml/reactions.md | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/yaml/reactions.md b/doc/sphinx/yaml/reactions.md index 93c98cff296..44496a4a5ae 100644 --- a/doc/sphinx/yaml/reactions.md +++ b/doc/sphinx/yaml/reactions.md @@ -118,7 +118,10 @@ a mapping with fields: : The pre-exponential factor `b` -: The temperature exponent, which is applied to the electron temperature +: The electron temperature exponent + +`b-gas` +: The gas temperature exponent `Ea-gas` : The activation energy term $E_{a,g}$ that is related to the gas temperature @@ -129,7 +132,7 @@ a mapping with fields: or a corresponding four-element list. The following are equivalent: ```yaml -{A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol} +{A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol, b-gas: 0} [17283, -3.1, -5820 J/mol, 1081 J/mol] ``` From 085147d2a2ea5fc2b726aaea84cf1fcf24b3cd5b Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 14:13:00 +0200 Subject: [PATCH 06/27] simplify twoTempPlasmaRate anymap constructor code --- src/kinetics/TwoTempPlasmaRate.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index 5f257fbb0a9..d97bc4cd88c 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -71,10 +71,7 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_u : TwoTempPlasmaRate() { setParameters(node, rate_units); - - if (node.hasKey("b-gas")) { - m_bg = node["b-gas"].asDouble(); - } + m_bg = node.getDouble("b-gas", 0.0); } double TwoTempPlasmaRate::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const From 6260a5a27695d4a023fd3d776bc9d032ee4a9cf2 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 15:00:57 +0200 Subject: [PATCH 07/27] updated cantera.bib --- doc/doxygen/cantera.bib | 59 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/doc/doxygen/cantera.bib b/doc/doxygen/cantera.bib index 6e9760f8fa1..d983b6b1de5 100644 --- a/doc/doxygen/cantera.bib +++ b/doc/doxygen/cantera.bib @@ -31,6 +31,22 @@ @article{blowers2000 url = {https://dx.doi.org/10.1002/aic.690461015}, volume = {46}, year = {2000}} +@book{capitelli2013, + title={Plasma kinetics in atmospheric gases}, + author={Capitelli, Mario and Ferreira, Carlos M and Gordiets, Boris F and Osipov, Alexey I}, + volume={31}, + year={2013}, + publisher={Springer Science \& Business Media} +} +@article{castela2016, + title={Modelling the impact of non-equilibrium discharges on reactive mixtures for simulations of plasma-assisted ignition in turbulent flows}, + author={Castela, Maria and Fiorina, Beno{\^\i}t and Coussement, Axel and Gicquel, Olivier and Darabiha, Nasser and Laux, Christophe O}, + journal={Combustion and flame}, + volume={166}, + pages={133--147}, + year={2016}, + publisher={Elsevier} +} @book{chapman1970, title={The mathematical theory of non-uniform gases}, author={S.~Chapman and T.~G.~Cowling}, @@ -111,6 +127,16 @@ @misc{gou2011 howpublished = {\url{https://engine.princeton.edu/model-reduction/}}, note = {Accessed: 2022-04-25}, year = {2011}} +@article{guerra2019, + title={Modelling N2--O2 plasmas: volume and surface kinetics}, + author={Guerra, Vasco and Tejero-del-Caz, Antonio and Pintassilgo, Carlos D and Alves, Lu{\'\i}s L}, + journal={Plasma Sources Science and Technology}, + volume={28}, + number={7}, + pages={073001}, + year={2019}, + publisher={IOP Publishing} +} @article{gudmundsson2001, author = {J.~T.~Gudmundsson}, journal = {Plasma Sources Science and Technology}, @@ -390,6 +416,16 @@ @book{poling2001 isbn = {0-07-149999-7}, year = {2001}, edition = {Fifth}} +@article{popov2016, + title={Kinetics of plasma-assisted combustion: effect of non-equilibrium excitation on the ignition and oxidation of combustible mixtures}, + author={Popov, NA}, + journal={Plasma Sources Science and Technology}, + volume={25}, + number={4}, + pages={043002}, + year={2016}, + publisher={IOP Publishing} +} @phdthesis{prager2005, author = {J.~Prager}, school = {Technische Universität Darmstadt}, @@ -470,6 +506,16 @@ @article{singal2024 number = {1--4}, pages = {105779}, year = {2024}} +@article{starikovskiy2013, +title = "Plasma-assisted ignition and combustion", +journal = "Progress in Energy and Combustion Science", +volume = "39", +number = "1", +pages = "61 - 110", +year = "2013", +issn = "0360-1285", +author = "Andrey Starikovskiy and Nickolay Aleksandrov", +} @article{stewart1989, author = {P.~H.~Stewart and C.~W.~Larson and D.~Golden}, journal = {Combustion and Flame}, @@ -570,6 +616,19 @@ @article{westbrook1981 doi = {10.1080/00102208108946970}, URL = {https://doi.org/10.1080/00102208108946970}, year = {1981}} +@article{zhong2023, + author = {Hongtao Zhong and Xingqian Mao and Ning Liu and Ziyu Wang and Timothy Ombrello and Yiguang Ju}, + doi = {10.1016/j.combustflame.2023.112948}, + issn = {15562921}, + pages = {112948}, + journal = {Combust. Flame}, + keywords = {Ammonia,Laser diagnostics,N2O/NOx chemistry,Nitrogen cycle,Plasmas,Sustainable energy}, + month = {10}, + publisher = {Elsevier Inc.}, + title = {{Understanding non-equilibrium N$_2$O/NO$_x$ chemistry in plasma-assisted low-temperature NH$_3$ oxidation}}, + volume = {256}, + year = {2023} +} @article{zhu2006, author = {H.~Zhu and R.~J.~Kee}, journal = {Journal of The Electrochemical Society}, From f31351dae77d485bf2c8c0449d1ec8d691d2d2e3 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 15:37:58 +0200 Subject: [PATCH 08/27] update documentation content and format in VibrationalRelaxationRate.h --- .../kinetics/VibrationalRelaxationRate.h | 141 ++++++++---------- 1 file changed, 62 insertions(+), 79 deletions(-) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index 2df4def6e8a..1d61e8d3af7 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -1,49 +1,10 @@ //! @file VibrationalRelaxationRate.h //! Header for vibrational relaxation reaction rates in plasma kinetics. -//! +//! @since New in %Cantera 4.0 + //! This file is part of Cantera. See License.txt in the top-level directory or //! at https://cantera.org/license.txt for license and copyright information. - -//! It implements the relaxation rate of vibrationally excited species in plasma kinetics. -//! Four models are currently supported by this class: 'constant', 'multi-state-resolved', 'starikovskiy', and 'castela'. - -//! The 'constant' model relaxes the vibrational species with a constant rate coefficient. It could just as well be an -//! Arrhenius rate, but the constant model is provided for convenience and to avoid confusion with conventional Arrhenius rates. - -//! The 'multi-state-resolved' model is meant to fully resolve vibrational relaxation by taking into account all vibrational species -//! in the phase (Ex: N2(v=1-8)) and solves for their V-T and V-V relaxation. The scaling of the rates are based on the SSH theory -//! detailed in Chapter 7 of the following book: -//! Capitelli, M., Ferreira, C. M., Gordiets, B. F., & Osipov, A. I. (2013). -//! Plasma kinetics in atmospheric gases (Vol. 31). Springer Science & Business Media. -//! The simplified version of the SSH theory implemented in this model is based on the harmonic oscillator approximation and can be found -//! in the following paper (equations 18 and 19): -//! Vasco Guerra et al 2019 Plasma Sources Sci. Technol. 28 073001 -//! The k10 rates are taken from //! Zhong, H. et al. (2023) “Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted -//! low-temperature NH3 oxidation,” Combustion and Flame, 256, which itself took those rates from both the Capitelli book and -//! the article from Strikovskiy and Aleksandrov (2013), see more below. - -//! The 'castela' model is meant to be used only for N2 vibrational relaxation, by collisions with N2, O2 and O exclusively. -//! It implements the mean vibrational energy relaxation model using a fictive cantera species. -//! The Castela model is based on the following paper: -//! Castela, M., Fiorina, B., Coussement, A., Gicquel, O., Darabiha, N., & Laux, C. O. (2016). -//! Modelling the impact of non-equilibrium discharges on reactive mixtures for simulations of -//! plasma-assisted ignition in turbulent flows. Combustion and flame, 166, 133-147. - -//! The so-called 'starikovskiy' model is an extension of the Castela model to several vibrational species and additional colliders. -//! A lot of vibrational relaxation rates can be found in Table 1 of the following paper: -//! Starikovskiy, A., & Aleksandrov, N. (2013). Plasma-assisted ignition and combustion. -//! Progress in Energy and Combustion Science, 39(1), 61-110. -//! The rates for the vibrational relaxation of NH3 can be found in the reaction mechanism provided in supplementary material -//! of the following paper: -//! Zhong, H. et al. (2023) “Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted low-temperature NH3 oxidation,” -//! Combustion and Flame, 256. -//! More rates for the vibrational relaxation of CH$ can be found in the following paper: -//! Popov, N.A. (2016) “Kinetics of plasma-assisted combustion: Effect of non-equilibrium excitation on the ignition -//! and oxidation of combustible mixtures,” Plasma Sources Science and Technology. Institute of Physics Publishing. - - - #ifndef CT_VIBRATIONALRELAXATIONRATE_H #define CT_VIBRATIONALRELAXATIONRATE_H @@ -53,12 +14,11 @@ namespace Cantera { - //! Shared temperature data for vibrational relaxation rates. struct DetailedVibData : public ReactionData { //! Update cached temperature-dependent data. - /*! + /** * @param phase Thermodynamic phase used to retrieve the gas temperature. * @param kin Kinetics object. Not used here, but required by the * ReactionData interface. @@ -71,9 +31,8 @@ struct DetailedVibData : public ReactionData using ReactionData::update; }; - -//! Vibrational relaxation reaction rate. -/*! +//! Vibrational relaxation reaction rate class definition. +/** * This class provides a common implementation for several vibrational * relaxation models: * @@ -107,12 +66,31 @@ struct DetailedVibData : public ReactionData * The selected physical model is specified separately using: * * @code{.yaml} - * vibration_model: constant + * type: vibrational-relaxation * vibration_model: multi-state-resolved - * vibration_model: starikovskiy - * vibration_model: castela * @endcode * + * Accepted values for `vibration_model` are `constant`, + * `multi-state-resolved`, `starikovskiy`, and `castela`. + * + * The `constant` model relaxes the vibrational species with a constant rate coefficient. It could just as well be an + * Arrhenius rate, but the constant model is provided for convenience and to avoid confusion with conventional Arrhenius rates + * for the YAML file user. + * + * The `multi-state-resolved` model is meant to fully resolve vibrational relaxation by taking into account all vibrational species + * in the phase (for example: N2(v=1-8)) and solves for their V-T and V-V relaxation. The scaling of the rates is based on the SSH theory + * detailed in Chapter 7 of @cite capitelli2013. The simplified version of the SSH theory implemented in this model is based on the + * harmonic oscillator approximation and can be found in equations 18 and 19 of @cite guerra2019. + * The @f$ k_{10} @f$ rates are taken from @cite zhong2023, @cite capitelli2013 and @cite starikovskiy2013. + * + * The `castela` model is meant to be used only for N2 vibrational relaxation, by collisions with N2, O2 and O exclusively. + * It implements the mean vibrational energy relaxation model using a fictitious Cantera species and is based on @cite castela2016. + * + * The model coined `starikovskiy` is an extension of the Castela model to several vibrational species and additional colliders. + * Many vibrational relaxation rates can be found in Table 1 of @cite starikovskiy2013 (hence the model name). + * The rates for the vibrational relaxation of NH3 can be found in the reaction mechanism provided in supplementary material of @cite zhong2023. + * More rates for the vibrational relaxation of @f$ \mathrm{CH_4} @f$ can be found in @cite popov2016. + * * Unit conventions: * * - `A` uses standard Cantera rate coefficient units. Its units depend on the @@ -125,7 +103,9 @@ struct DetailedVibData : public ReactionData * The coefficients `B`, `C`, `D`, `E`, `m`, `z`, and `scaling` are read as * raw floating-point values. They are not converted by Cantera's unit system. * + * * @ingroup arrheniusGroup + * @since New in %Cantera 4.0 */ class VibrationalRelaxationRate : public ArrheniusBase { @@ -134,7 +114,7 @@ class VibrationalRelaxationRate : public ArrheniusBase VibrationalRelaxationRate(); //! Constructor using the internal generic representation. - /*! + /** * @param A Pre-exponential factor. * @param B Dimensionless constant in the exponential. * @param C Coefficient multiplying T^(-1/3). @@ -144,43 +124,49 @@ class VibrationalRelaxationRate : public ArrheniusBase * @param m Temperature exponent used by the D term. * @param E Coefficient multiplying T^(-z). * @param z Temperature exponent used by the E term. + * @since New in %Cantera 4.0 */ VibrationalRelaxationRate(double A, double B, double C, double D, - double b, double scaling = 1.0, - double m = 2.0 / 3.0, - double E = 0.0, double z = 1.0); + double b, double scaling = 1.0, + double m = 2.0 / 3.0, + double E = 0.0, double z = 1.0); //! Constructor based on AnyMap content. - explicit VibrationalRelaxationRate(const AnyMap& node, - const UnitStack& rate_units = {}); + //! @since New in %Cantera 4.0 + explicit VibrationalRelaxationRate(const AnyMap& node, const UnitStack& rate_units = {}); //! Set rate parameters from an AnyMap. - void setParameters(const AnyMap& node, - const UnitStack& rate_units) override; + //! @since New in %Cantera 4.0 + void setParameters(const AnyMap& node, const UnitStack& rate_units) override; //! Get rate parameters for YAML serialization. + //! @since New in %Cantera 4.0 void getParameters(AnyMap& node) const override; //! Create a rate evaluator for this reaction rate type. + //! @since New in %Cantera 4.0 unique_ptr newMultiRate() const override { return make_unique>(); } //! String identifying this reaction rate type. + //! @since New in %Cantera 4.0 const string type() const override { return "vibrational-relaxation"; } //! Set context of reaction rate evaluation. - /*! + /** * Vibrational relaxation rates are intended for irreversible * non-equilibrium plasma reactions. Reversible reactions are rejected * because the reverse rate cannot be obtained from conventional * thermochemistry for these models. */ + //! @since New in %Cantera 4.0 void setContext(const Reaction& rxn, const Kinetics& kin) override; //! Evaluate the forward rate coefficient. + //! @since New in %Cantera 4.0 double evalFromStruct(const DetailedVibData& shared_data) const { const double invT = shared_data.recipT; const double invT13 = std::cbrt(invT); @@ -195,7 +181,7 @@ class VibrationalRelaxationRate : public ArrheniusBase } //! Evaluate the scaled temperature derivative. - /*! + /** * This returns: * * @f[ @@ -213,6 +199,7 @@ class VibrationalRelaxationRate : public ArrheniusBase * - z E T^{-z-1} * @f] */ + //! @since New in %Cantera 4.0 double ddTScaledFromStruct(const DetailedVibData& shared_data) const; private: @@ -247,9 +234,8 @@ class VibrationalRelaxationRate : public ArrheniusBase double m_referencePressure = OneAtm; //! Selected vibrational relaxation model. - /*! + /** * Accepted values: - * * - `constant` * - `multi-state-resolved` * - `starikovskiy` @@ -257,7 +243,7 @@ class VibrationalRelaxationRate : public ArrheniusBase */ string m_vibration_model = "multi-state-resolved"; - //! YAML key names. + //! YAML variable names. string m_B_str = "B"; string m_C_str = "C"; string m_D_str = "D"; @@ -268,26 +254,23 @@ class VibrationalRelaxationRate : public ArrheniusBase string m_reference_pressure_str = "reference-pressure"; //! Configure the ArrheniusBase part from an already-converted internal A value. - /*! - * This is needed for models such as Castela, where the user-facing YAML does - * not contain a standard Arrhenius A coefficient. - */ - void configureBaseFromInternalA(const AnyMap& node, - const UnitStack& rate_units, + /** + * This is needed for models such as Castela, where the user-facing YAML does + * not contain a standard Arrhenius A coefficient. + * @since New in %Cantera 4.0 + */ + void configureBaseFromInternalA(const AnyMap& node, const UnitStack& rate_units, double A, double b); //! Configure the ArrheniusBase part from a YAML A value and an explicit b. - /*! - * This is needed for models such as constant and Starikovskiy, where the YAML - * does not contain the standard Arrhenius pair A / b. - */ - void configureBaseFromYamlA(const AnyMap& node, - const UnitStack& rate_units, - const AnyValue& A, - double b); - - }; + /** + * This is needed for models such as `constant` and `starikovskiy`, where the YAML + * does not contain the standard Arrhenius pair A / b. + * @since New in %Cantera 4.0 + */ + void configureBaseFromYamlA(const AnyMap& node, const UnitStack& rate_units, + const AnyValue& A, double b); +}; } - #endif \ No newline at end of file From ad4f41d3651c1e5a0d27d6ca68bc0efbb0dfc3d2 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 15:46:51 +0200 Subject: [PATCH 09/27] line length formatting + further formatting of VibrationalRelaxationRate.h --- .../kinetics/VibrationalRelaxationRate.h | 71 ++++++++++--------- 1 file changed, 36 insertions(+), 35 deletions(-) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index 1d61e8d3af7..95ebc93c4e4 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -1,7 +1,7 @@ //! @file VibrationalRelaxationRate.h //! Header for vibrational relaxation reaction rates in plasma kinetics. //! @since New in %Cantera 4.0 - +//! //! This file is part of Cantera. See License.txt in the top-level directory or //! at https://cantera.org/license.txt for license and copyright information. @@ -22,7 +22,6 @@ struct DetailedVibData : public ReactionData * @param phase Thermodynamic phase used to retrieve the gas temperature. * @param kin Kinetics object. Not used here, but required by the * ReactionData interface. - * * @return `true` if the temperature has changed and rates need to be * recomputed; `false` otherwise. */ @@ -35,14 +34,11 @@ struct DetailedVibData : public ReactionData /** * This class provides a common implementation for several vibrational * relaxation models: - * * - `constant` * - `multi-state-resolved` * - `starikovskiy` * - `castela` - * * Internally, all models are mapped to the following generic expression: - * * @f[ * k_f = * scaling \, A \, @@ -54,43 +50,47 @@ struct DetailedVibData : public ReactionData * + E T^{-z} * \right) * @f] - * * where `T` is the gas temperature in K. * * The YAML reaction type is: - * * @code{.yaml} * type: vibrational-relaxation * @endcode - * * The selected physical model is specified separately using: - * * @code{.yaml} * type: vibrational-relaxation * vibration_model: multi-state-resolved * @endcode - * * Accepted values for `vibration_model` are `constant`, * `multi-state-resolved`, `starikovskiy`, and `castela`. - * - * The `constant` model relaxes the vibrational species with a constant rate coefficient. It could just as well be an - * Arrhenius rate, but the constant model is provided for convenience and to avoid confusion with conventional Arrhenius rates - * for the YAML file user. - * - * The `multi-state-resolved` model is meant to fully resolve vibrational relaxation by taking into account all vibrational species - * in the phase (for example: N2(v=1-8)) and solves for their V-T and V-V relaxation. The scaling of the rates is based on the SSH theory - * detailed in Chapter 7 of @cite capitelli2013. The simplified version of the SSH theory implemented in this model is based on the - * harmonic oscillator approximation and can be found in equations 18 and 19 of @cite guerra2019. - * The @f$ k_{10} @f$ rates are taken from @cite zhong2023, @cite capitelli2013 and @cite starikovskiy2013. - * - * The `castela` model is meant to be used only for N2 vibrational relaxation, by collisions with N2, O2 and O exclusively. - * It implements the mean vibrational energy relaxation model using a fictitious Cantera species and is based on @cite castela2016. * - * The model coined `starikovskiy` is an extension of the Castela model to several vibrational species and additional colliders. - * Many vibrational relaxation rates can be found in Table 1 of @cite starikovskiy2013 (hence the model name). - * The rates for the vibrational relaxation of NH3 can be found in the reaction mechanism provided in supplementary material of @cite zhong2023. - * More rates for the vibrational relaxation of @f$ \mathrm{CH_4} @f$ can be found in @cite popov2016. - * + * The `constant` model relaxes the vibrational species with a constant rate + * coefficient. It could just as well be an Arrhenius rate, but the constant + * model is provided for convenience and to avoid confusion with conventional + * Arrhenius rates for the YAML file user. + * + * The `multi-state-resolved` model fully resolves vibrational relaxation by + * taking into account all vibrational species in the phase, for example + * `N2(v=1-8)`, and solves for their V-T and V-V relaxation. The scaling of the + * rates is based on the SSH theory detailed in Chapter 7 of + * @cite capitelli2013. The simplified SSH theory implemented here is based on + * the harmonic oscillator approximation and can be found in equations 18 and + * 19 of @cite guerra2019. The @f$ k_{10} @f$ rates are taken from + * @cite zhong2023, @cite capitelli2013, and @cite starikovskiy2013. + * + * The `castela` model is meant to be used only for N2 vibrational relaxation, + * by collisions with N2, O2, and O exclusively. It implements the mean + * vibrational energy relaxation model using a fictitious Cantera species and + * is based on @cite castela2016. + * + * The `starikovskiy` model is an extension of the Castela model to several + * vibrational species and additional colliders. Many vibrational relaxation + * rates can be found in Table 1 of @cite starikovskiy2013, hence the model + * name. The rates for the vibrational relaxation of NH3 can be found in the + * reaction mechanism provided in the supplementary material of + * @cite zhong2023. More rates for the vibrational relaxation of + * @f$ \mathrm{CH_4} @f$ can be found in @cite popov2016. + * * Unit conventions: * * - `A` uses standard Cantera rate coefficient units. Its units depend on the @@ -103,7 +103,6 @@ struct DetailedVibData : public ReactionData * The coefficients `B`, `C`, `D`, `E`, `m`, `z`, and `scaling` are read as * raw floating-point values. They are not converted by Cantera's unit system. * - * * @ingroup arrheniusGroup * @since New in %Cantera 4.0 */ @@ -133,7 +132,8 @@ class VibrationalRelaxationRate : public ArrheniusBase //! Constructor based on AnyMap content. //! @since New in %Cantera 4.0 - explicit VibrationalRelaxationRate(const AnyMap& node, const UnitStack& rate_units = {}); + explicit VibrationalRelaxationRate(const AnyMap& node, + const UnitStack& rate_units = {}); //! Set rate parameters from an AnyMap. //! @since New in %Cantera 4.0 @@ -146,7 +146,8 @@ class VibrationalRelaxationRate : public ArrheniusBase //! Create a rate evaluator for this reaction rate type. //! @since New in %Cantera 4.0 unique_ptr newMultiRate() const override { - return make_unique>(); + return make_unique>(); } //! String identifying this reaction rate type. @@ -259,16 +260,16 @@ class VibrationalRelaxationRate : public ArrheniusBase * not contain a standard Arrhenius A coefficient. * @since New in %Cantera 4.0 */ - void configureBaseFromInternalA(const AnyMap& node, const UnitStack& rate_units, + void configureBaseFromInternalA(const AnyMap& node, const UnitStack& rate_units, double A, double b); //! Configure the ArrheniusBase part from a YAML A value and an explicit b. /** - * This is needed for models such as `constant` and `starikovskiy`, where the YAML - * does not contain the standard Arrhenius pair A / b. + * This is needed for models such as `constant` and `starikovskiy`, + * where the YAML does not contain the standard Arrhenius pair A / b. * @since New in %Cantera 4.0 */ - void configureBaseFromYamlA(const AnyMap& node, const UnitStack& rate_units, + void configureBaseFromYamlA(const AnyMap& node, const UnitStack& rate_units, const AnyValue& A, double b); }; From 54d996be6a81f7fcfabb1be179922374e095ae1e Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 16:12:40 +0200 Subject: [PATCH 10/27] modify basic formatting in vibrationalRelaxationRate.cpp --- src/kinetics/VibrationalRelaxationRate.cpp | 30 +++------------------- 1 file changed, 3 insertions(+), 27 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 79aabf911ee..6a11aa830c4 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -1,4 +1,5 @@ //! @file VibrationalRelaxationRate.cpp +//! @since New in %Cantera 4.0 // This file is part of Cantera. See License.txt in the top-level directory or // at https://cantera.org/license.txt for license and copyright information. @@ -20,7 +21,6 @@ namespace Cantera namespace { - void requireNoKey(const AnyMap& node, const string& key, const string& model, const string& where) { @@ -30,7 +30,6 @@ void requireNoKey(const AnyMap& node, const string& key, } } - void requireKey(const AnyMap& node, const string& key, const string& model, const string& where) { @@ -40,22 +39,20 @@ void requireKey(const AnyMap& node, const string& key, } } - bool isVibrationalSpecies(const string& name) { - // Supported names: + // Supported names formats: // // N2(v) // O2(v0) // O2(v1) // O2(v12) // - // The old form O2(v=1) is intentionally not required. + // Style like O2(v=1) should not be used. const auto pos = name.find("(v"); return pos != string::npos && !name.empty() && name.back() == ')'; } - string groundStateName(const string& name) { const auto pos = name.find("(v"); @@ -65,14 +62,12 @@ string groundStateName(const string& name) return name.substr(0, pos); } - string vibrationalFamilyName(const string& name) { // Collapse O2(v1), O2(v2), O2(v12), and O2(v) into O2(v). return groundStateName(name) + "(v)"; } - double compositionSum(const Composition& comp) { double sum = 0.0; @@ -82,7 +77,6 @@ double compositionSum(const Composition& comp) return sum; } - Composition replaceVibrationalSpeciesByGroundState(const Composition& comp) { Composition out; @@ -99,7 +93,6 @@ Composition replaceVibrationalSpeciesByGroundState(const Composition& comp) return out; } - bool sameComposition(const Composition& a, const Composition& b, double tol = 1e-12) { @@ -134,7 +127,6 @@ bool sameComposition(const Composition& a, const Composition& b, return true; } - std::vector vibrationalSpeciesInComposition(const Composition& comp) { std::vector out; @@ -151,7 +143,6 @@ std::vector vibrationalSpeciesInComposition(const Composition& comp) return out; } - // Registry used to check that a given vibrational family uses exactly one // relaxation model inside one Kinetics object. // @@ -169,7 +160,6 @@ std::vector vibrationalSpeciesInComposition(const Composition& comp) std::map> s_modelByFamily; std::set s_warnedMixedModels; - void registerVibrationalModelConsistency(const Kinetics& kin, const string& family, const string& model, @@ -209,7 +199,6 @@ void registerVibrationalModelConsistency(const Kinetics& kin, } } - string inferRelaxingFamily(const Reaction& rxn) { const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); @@ -223,7 +212,6 @@ string inferRelaxingFamily(const Reaction& rxn) return vibrationalFamilyName(vibReactants.front()); } - void validateSimpleRelaxationToGroundState(const Reaction& rxn, const string& model) { @@ -261,7 +249,6 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, } } - void validateCastelaReaction(const Reaction& rxn) { validateSimpleRelaxationToGroundState(rxn, "castela"); @@ -301,7 +288,6 @@ void validateCastelaReaction(const Reaction& rxn) } } - void validateDetailedRelaxationReaction(const Reaction& rxn) { const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); @@ -365,7 +351,6 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) } // namespace - bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) { const double T = phase.temperature(); @@ -378,12 +363,10 @@ bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) return true; } - VibrationalRelaxationRate::VibrationalRelaxationRate() { } - VibrationalRelaxationRate::VibrationalRelaxationRate( double A, double B, double C, double D, double b, double scaling, double m, double E, double z) @@ -398,7 +381,6 @@ VibrationalRelaxationRate::VibrationalRelaxationRate( { } - VibrationalRelaxationRate::VibrationalRelaxationRate( const AnyMap& node, const UnitStack& rate_units) : VibrationalRelaxationRate() @@ -406,7 +388,6 @@ VibrationalRelaxationRate::VibrationalRelaxationRate( setParameters(node, rate_units); } - void VibrationalRelaxationRate::configureBaseFromInternalA( const AnyMap& node, const UnitStack& rate_units, double A, double b) { @@ -433,7 +414,6 @@ void VibrationalRelaxationRate::configureBaseFromInternalA( m_valid = true; } - void VibrationalRelaxationRate::configureBaseFromYamlA( const AnyMap& node, const UnitStack& rate_units, const AnyValue& A, double b) @@ -460,7 +440,6 @@ void VibrationalRelaxationRate::configureBaseFromYamlA( m_valid = true; } - void VibrationalRelaxationRate::setParameters(const AnyMap& node, const UnitStack& rate_units) { @@ -696,7 +675,6 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, } } - void VibrationalRelaxationRate::getParameters(AnyMap& node) const { if (!valid()) { @@ -791,7 +769,6 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const node["rate-constant"] = std::move(rateNode); } - double VibrationalRelaxationRate::ddTScaledFromStruct( const DetailedVibData& shared_data) const { @@ -804,7 +781,6 @@ double VibrationalRelaxationRate::ddTScaledFromStruct( - m_z * m_E * std::pow(invT, m_z) * invT; } - void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& kin) { if (rxn.reversible) { From 1165179f968351139364274c9b4ae108175f9e01 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 18:00:30 +0200 Subject: [PATCH 11/27] Code formatting and commenting in Plasma phase .h and .cpp --- include/cantera/thermo/PlasmaPhase.h | 17 +++++++++-- src/thermo/PlasmaPhase.cpp | 43 ++++++++++++++++++---------- 2 files changed, 43 insertions(+), 17 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 0b3b1e2bce4..f0803b4ebaa 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -1026,19 +1026,32 @@ class PlasmaPhase: public IdealGasPhase //! Work array mutable std::vector m_work; + //! A structure to describe species serving as mean vibrational energy reservoirs struct VibrationalReservoirSpecies { - size_t reservoirIndex = npos; - size_t baseSpeciesIndex = npos; + size_t reservoirIndex = npos; // index of the vibrational reservoir species + size_t baseSpeciesIndex = npos; // index of the corresponding ground state phase species }; + //! Vector of species serving as mean vibrational energy reservoirs std::vector m_vibrationalReservoirSpecies; + //! A boolean flag to update vibrational reservoir species bool m_vibrationalReservoirSpeciesNeedUpdate = true; + //! Threshold fraction of vibrational reservoirs to their ground state + //! above which they may be at risk of hindering phase chemistry. double m_vibrationalMoleFractionThreshold = 1e-2; + + //! The absolute mole threshold below which a the chemistry is + //! assumed to be safe from alterations from vibrational reservoirs. double m_vibrationalAbsoluteMoleFractionThreshold = 1e-20; + //! Detects species serving as vibrational reservoirs and + //! adds them to m_vibrationalReservoirSpecies. void updateVibrationalReservoirSpecies(); + + //! A function to check that vibrational reservoir species + //! are not at risk to hinder phase chemistry. void checkVibrationalReservoirMoleFractions(); }; diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 69c6cab2cfd..deedd3f8afe 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -20,6 +20,9 @@ namespace Cantera { namespace { const double gamma = sqrt(2 * ElectronCharge / ElectronMass); + //! A function to check whether a phase species is actually a vibrational reservoir species + //! If the species is actually a reservoir species, the function sets the basename + //! at the given baseName adress. bool isVibrationalReservoirName(const string& name, string& baseName) { const string suffix = "(v)"; @@ -126,6 +129,8 @@ bool PlasmaPhase::addSpecies(shared_ptr spec) } } + // Adding species may introduce new `X(v)` reservoirs or their base species. + // Set the flag to true to run a check. if (added) { m_vibrationalReservoirSpeciesNeedUpdate = true; } @@ -1027,6 +1032,7 @@ double PlasmaPhase::intrinsicHeating() const double qJ = jouleHeatingPower(); checkFinite(qJ); + // set the check here to be updated at runtime checkVibrationalReservoirMoleFractions(); return qJ; @@ -1046,6 +1052,7 @@ void PlasmaPhase::updateVibrationalReservoirSpecies() size_t kBase = speciesIndex(baseName, false); + // Check that the base species of the vibrational reservoir actually exists in the phase. if (kBase == npos) { warn_user("PlasmaPhase::updateVibrationalReservoirSpecies", "Species '{}' matches the fictive vibrational-reservoir naming " @@ -1061,7 +1068,7 @@ void PlasmaPhase::updateVibrationalReservoirSpecies() m_vibrationalReservoirSpecies.push_back(reservoir); } - + // the update being done, the flag can be set back to false. m_vibrationalReservoirSpeciesNeedUpdate = false; } @@ -1075,33 +1082,39 @@ void PlasmaPhase::checkVibrationalReservoirMoleFractions() const size_t kReservoir = reservoir.reservoirIndex; const size_t kBase = reservoir.baseSpeciesIndex; - const double Xv = moleFraction(kReservoir); - const double Xb = moleFraction(kBase); + const double Xv = moleFraction(kReservoir); // mole fraction of the vibrational reservoir species + const double Xb = moleFraction(kBase); // mole fraction of the base species it is associated to const double pool = Xv + Xb; + // Ignore species pools that are too diluted to meaningfully affect chemistry. if (pool <= m_vibrationalAbsoluteMoleFractionThreshold) { continue; } + // Check that the fraction of vibrational species is not too high + // with respect to its base species. Should this fraction be too high, there + // is a risk for the phase chemistry to be altered by the reservoir: + // the code raises a warning to the user. const double reservoirFraction = Xv / pool; if (reservoirFraction > m_vibrationalMoleFractionThreshold) { const string& reservoirName = speciesName(kReservoir); const string& baseName = speciesName(kBase); - writelog("Warning: fictive vibrational reservoir species '{}' contains " - "{:.3e} of the total '{}' pool. " - "X({}) = {:.3e}, X({}) = {:.3e}, threshold = {:.3e}. " - "Chemistry involving '{}' may be affected because part of the " - "material is stored in an inert vibrational reservoir.\n", - reservoirName, - reservoirFraction, - baseName, - reservoirName, Xv, - baseName, Xb, - m_vibrationalMoleFractionThreshold, - baseName); + warn_user("PlasmaPhase::checkVibrationalReservoirMoleFractions", + "Warning: fictive vibrational reservoir species '{}' contains " + "{:.3e} of the total '{}' pool. " + "X({}) = {:.3e}, X({}) = {:.3e}, threshold = {:.3e}. " + "Chemistry involving '{}' may be affected because part of the " + "material is stored in an inert vibrational reservoir.\n", + reservoirName, + reservoirFraction, + baseName, + reservoirName, Xv, + baseName, Xb, + m_vibrationalMoleFractionThreshold, + baseName); } } } From dc3580dda6df29d376985243dbc4561119242224 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 19:51:29 +0200 Subject: [PATCH 12/27] Add comments and simplify some code bits in vibrationalRelaxationRate.cpp --- src/kinetics/VibrationalRelaxationRate.cpp | 73 +++++++++++----------- 1 file changed, 38 insertions(+), 35 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 6a11aa830c4..bdd6d6b4a31 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -21,6 +21,7 @@ namespace Cantera namespace { +// helpers to check the correctness of chosen vibrational model. void requireNoKey(const AnyMap& node, const string& key, const string& model, const string& where) { @@ -39,6 +40,8 @@ void requireKey(const AnyMap& node, const string& key, } } +// check whether a species is a vibrational excited species +// or a vibrational reservoir from its name bool isVibrationalSpecies(const string& name) { // Supported names formats: @@ -53,6 +56,7 @@ bool isVibrationalSpecies(const string& name) return pos != string::npos && !name.empty() && name.back() == ')'; } +// find the name of the ground state of a vibrationally excited species. string groundStateName(const string& name) { const auto pos = name.find("(v"); @@ -64,7 +68,8 @@ string groundStateName(const string& name) string vibrationalFamilyName(const string& name) { - // Collapse O2(v1), O2(v2), O2(v12), and O2(v) into O2(v). + // Collapse O2(v1), O2(v2), O2(v12), and O2(v) into O2(v) as a "family + // name" to identify which species have the same ground state "parent". return groundStateName(name) + "(v)"; } @@ -146,13 +151,15 @@ std::vector vibrationalSpeciesInComposition(const Composition& comp) // Registry used to check that a given vibrational family uses exactly one // relaxation model inside one Kinetics object. // -// Allowed: +// It is allowed for each "vibrational family" +// to have its own vibrational relaxation model: // // N2(v) -> constant // O2(v) -> multi-state-resolved // NH3(v) -> starikovskiy // -// Forbidden: +// But it is forbidden for each reaction of the same "vibrational family" +// to have different relaxation models: // // N2(v) + O -> castela // N2(v) + N2 -> starikovskiy @@ -300,7 +307,7 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) } // All vibrational species in a detailed relaxation reaction are required - // to belong to the same vibrational family. + // to belong to the same "vibrational family". const string family = vibrationalFamilyName(vibReactants.front()); for (const auto& sp : vibReactants) { @@ -349,7 +356,7 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) } } -} // namespace +} // end of namespace where all the helpers and safety functions are defined. bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) { @@ -363,10 +370,12 @@ bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) return true; } +// Default constructor VibrationalRelaxationRate::VibrationalRelaxationRate() { } +// Constructor VibrationalRelaxationRate::VibrationalRelaxationRate( double A, double B, double C, double D, double b, double scaling, double m, double E, double z) @@ -381,6 +390,7 @@ VibrationalRelaxationRate::VibrationalRelaxationRate( { } +// AnyMap constructor VibrationalRelaxationRate::VibrationalRelaxationRate( const AnyMap& node, const UnitStack& rate_units) : VibrationalRelaxationRate() @@ -396,6 +406,7 @@ void VibrationalRelaxationRate::configureBaseFromInternalA( // We intentionally do not call ArrheniusBase::setParameters here because // some vibration models do not expose a standard YAML rate-constant with // both A and b. Castela is the main example. + // This is also done for future class extension compatibility. ReactionRate::setParameters(node, rate_units); setRateUnits(rate_units); @@ -410,7 +421,7 @@ void VibrationalRelaxationRate::configureBaseFromInternalA( m_logA = NAN; } - // Critical: ArrheniusBase::validate checks this flag. + // ArrheniusBase::validate checks this flag. m_valid = true; } @@ -436,7 +447,7 @@ void VibrationalRelaxationRate::configureBaseFromYamlA( m_logA = NAN; } - // Critical: ArrheniusBase::validate checks this flag. + // ArrheniusBase::validate checks this flag. m_valid = true; } @@ -449,19 +460,18 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // // The physical model is selected separately by: // - // vibration_model: constant - // vibration_model: multi-state-resolved - // vibration_model: starikovskiy - // vibration_model: castela + // vibration-model: constant + // vibration-model: multi-state-resolved + // vibration-model: starikovskiy + // vibration-model: castela // - // For backward compatibility, the default model is multi-state-resolved. - m_vibration_model = "multi-state-resolved"; + // This is done in a spriti of class adaptability: + // Any user wishing to develop more advanced vibration models + // is welcome to add them in this class. + // + // The default model is multi-state-resolved. - if (node.hasKey("vibration_model")) { - m_vibration_model = node["vibration_model"].asString(); - } else if (node.hasKey("model")) { - m_vibration_model = node["model"].asString(); - } + node.getString("vibration-model", "multi-state-resolved"); if (!node.hasKey("rate-constant")) { throw InputFileError("VibrationalRelaxationRate::setParameters", node, @@ -482,10 +492,8 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // Constant model: // // k(T) = A - requireKey(rateMap, m_A_str, m_vibration_model, "VibrationalRelaxationRate::setParameters"); - requireNoKey(rateMap, m_b_str, m_vibration_model, "VibrationalRelaxationRate::setParameters"); requireNoKey(rateMap, "n", m_vibration_model, @@ -527,7 +535,6 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, requireKey(rateMap, m_A_str, m_vibration_model, "VibrationalRelaxationRate::setParameters"); - requireNoKey(rateMap, "n", m_vibration_model, "VibrationalRelaxationRate::setParameters"); requireNoKey(rateMap, m_m_str, m_vibration_model, @@ -539,15 +546,13 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, ArrheniusBase::setParameters(node, rate_units); - m_B = rateMap.hasKey(m_B_str) ? rateMap[m_B_str].asDouble() : 0.0; - m_C = rateMap.hasKey(m_C_str) ? rateMap[m_C_str].asDouble() : 0.0; - m_D = rateMap.hasKey(m_D_str) ? rateMap[m_D_str].asDouble() : 0.0; + m_B = rateMap.getDouble(m_B_str, 0.0); + m_C = rateMap.getDouble(m_C_str, 0.0); + m_D = rateMap.getDouble(m_D_str, 0.0); m_m = 2.0 / 3.0; m_E = 0.0; m_z = 1.0; - m_scaling = rateMap.hasKey(m_scaling_str) - ? rateMap[m_scaling_str].asDouble() - : 1.0; + m_scaling = rateMap.getDouble(m_scaling_str, 1.0); } else if (m_vibration_model == "starikovskiy") { // User-facing formula: @@ -571,7 +576,6 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, requireKey(rateMap, m_A_str, m_vibration_model, "VibrationalRelaxationRate::setParameters"); - requireNoKey(rateMap, m_b_str, m_vibration_model, "VibrationalRelaxationRate::setParameters"); requireNoKey(rateMap, m_scaling_str, m_vibration_model, @@ -580,12 +584,12 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, const double n = rateMap.hasKey("n") ? rateMap["n"].asDouble() : 0.0; configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], n); - m_B = rateMap.hasKey("K") ? rateMap["K"].asDouble() : 0.0; - m_C = rateMap.hasKey("B") ? -rateMap["B"].asDouble() : 0.0; - m_D = rateMap.hasKey("C") ? rateMap["C"].asDouble() : 0.0; - m_m = rateMap.hasKey("m") ? rateMap["m"].asDouble() : 1.0; - m_E = rateMap.hasKey("D") ? rateMap["D"].asDouble() : 0.0; - m_z = rateMap.hasKey("z") ? rateMap["z"].asDouble() : 1.0; + m_B = rateMap.getDouble("K", 0.0); + m_C = rateMap.getDouble("B", 0.0); + m_D = rateMap.getDouble("C", 0.0); + m_m = rateMap.getDouble("m", 1.0); + m_E = rateMap.getDouble("D", 0.0); + m_z = rateMap.getDouble("z", 1.0); m_scaling = 1.0; if (m_m <= 0.0 || m_z <= 0.0) { @@ -620,7 +624,6 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, "VibrationalRelaxationRate::setParameters"); requireKey(rateMap, "b", m_vibration_model, "VibrationalRelaxationRate::setParameters"); - requireNoKey(rateMap, m_A_str, m_vibration_model, "VibrationalRelaxationRate::setParameters"); requireNoKey(rateMap, "n", m_vibration_model, From 0437f84182adadf5df71cfe037e2e9c8f1227bcb Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 20:02:03 +0200 Subject: [PATCH 13/27] correct mistake made on vibration-model assignment during previous commit --- src/kinetics/VibrationalRelaxationRate.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index bdd6d6b4a31..65654dcb409 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -471,7 +471,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // // The default model is multi-state-resolved. - node.getString("vibration-model", "multi-state-resolved"); + m_vibration_model = node.getString("vibration-model", "multi-state-resolved"); if (!node.hasKey("rate-constant")) { throw InputFileError("VibrationalRelaxationRate::setParameters", node, From c639bb20af35262094fe557e7dbe33ec82a6ff5c Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 22:06:31 +0200 Subject: [PATCH 14/27] drop starikovskiy model in-house sign convention to generalise the expression completely. The minus sign will be born by YAML coefficients directly. --- src/kinetics/VibrationalRelaxationRate.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 65654dcb409..7a63c4008b9 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -559,7 +559,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // // k(T) = A * T^n * exp( // K - // - B * T^(-1/3) + // + B * T^(-1/3) // + C * T^(-m) // + D * T^(-z) // ) @@ -568,7 +568,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // // m_b = n // m_B = K - // m_C = -B + // m_C = B // m_D = C // m_m = m // m_E = D @@ -731,7 +731,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const rateNode["n"] = m_b; rateNode["K"] = m_B; - rateNode["B"] = -m_C; + rateNode["B"] = m_C; rateNode["C"] = m_D; rateNode["m"] = m_m; rateNode["D"] = m_E; From 53a1605a431ca8e3271ad54fd0bf41c298b49c6f Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 22:20:05 +0200 Subject: [PATCH 15/27] replace all YAML parameters' underscores with hyphens --- .../kinetics/VibrationalRelaxationRate.h | 4 ++-- src/kinetics/VibrationalRelaxationRate.cpp | 20 +++++++++---------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index 95ebc93c4e4..825dad50a6a 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -59,9 +59,9 @@ struct DetailedVibData : public ReactionData * The selected physical model is specified separately using: * @code{.yaml} * type: vibrational-relaxation - * vibration_model: multi-state-resolved + * vibration-model: multi-state-resolved * @endcode - * Accepted values for `vibration_model` are `constant`, + * Accepted values for `vibration-model` are `constant`, * `multi-state-resolved`, `starikovskiy`, and `castela`. * * The `constant` model relaxes the vibrational species with a constant rate diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 7a63c4008b9..77c20cce392 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -27,7 +27,7 @@ void requireNoKey(const AnyMap& node, const string& key, { if (node.hasKey(key)) { throw InputFileError(where, node, - "Key '{}' is not allowed for vibration_model '{}'.", key, model); + "Key '{}' is not allowed for vibration-model '{}'.", key, model); } } @@ -36,7 +36,7 @@ void requireKey(const AnyMap& node, const string& key, { if (!node.hasKey(key)) { throw InputFileError(where, node, - "Missing required key '{}' for vibration_model '{}'.", key, model); + "Missing required key '{}' for vibration-model '{}'.", key, model); } } @@ -177,7 +177,7 @@ void registerVibrationalModelConsistency(const Kinetics& kin, const auto existing = modelByFamily.find(family); if (existing != modelByFamily.end() && existing->second != model) { throw InputFileError("VibrationalRelaxationRate::setContext", input, - "Inconsistent vibration_model for vibrational family '{}'. " + "Inconsistent vibration-model for vibrational family '{}'. " "This family was already registered with model '{}', but the " "current reaction uses model '{}'. A given vibrational family " "must use exactly one relaxation model.", @@ -227,13 +227,13 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, if (vibReactants.size() != 1) { throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, - "vibration_model '{}' expects exactly one vibrational reactant.", + "vibration-model '{}' expects exactly one vibrational reactant.", model); } if (!vibProducts.empty()) { throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, - "vibration_model '{}' describes relaxation to the ground state " + "vibration-model '{}' describes relaxation to the ground state " "and therefore does not allow vibrational products.", model); } @@ -251,7 +251,7 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) { throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, - "vibration_model '{}' expects a bimolecular relaxation reaction " + "vibration-model '{}' expects a bimolecular relaxation reaction " "of the form X(v) + M => X + M.", model); } } @@ -672,7 +672,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, } else { throw InputFileError("VibrationalRelaxationRate::setParameters", node, - "Unrecognized vibration_model '{}'. Expected 'multi-state-resolved', " + "Unrecognized vibration-model '{}'. Expected 'multi-state-resolved', " "'starikovskiy', 'castela', or 'constant'.", m_vibration_model); } @@ -688,7 +688,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const node["negative-A"] = true; } - node["vibration_model"] = m_vibration_model; + node["vibration-model"] = m_vibration_model; AnyMap rateNode; @@ -763,7 +763,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const } else { throw InputFileError("VibrationalRelaxationRate::getParameters", node, - "Unrecognized vibration_model '{}'. Expected 'multi-state-resolved', " + "Unrecognized vibration-model '{}'. Expected 'multi-state-resolved', " "'starikovskiy', 'castela', or 'constant'.", m_vibration_model); } @@ -804,7 +804,7 @@ void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& validateDetailedRelaxationReaction(rxn); } else { throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, - "Unrecognized vibration_model '{}'.", m_vibration_model); + "Unrecognized vibration-model '{}'.", m_vibration_model); } registerVibrationalModelConsistency( From 07807c98c1a391175a1cb4b5aab68f826caaeafa Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 22:31:34 +0200 Subject: [PATCH 16/27] centralise string handling for often used strings to avoid typos on them --- src/kinetics/VibrationalRelaxationRate.cpp | 157 +++++++++++---------- 1 file changed, 85 insertions(+), 72 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 77c20cce392..3243295a325 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -21,7 +21,20 @@ namespace Cantera namespace { -// helpers to check the correctness of chosen vibrational model. + +constexpr double Tiny = 1e-12; + +const string WhereSetParameters = "VibrationalRelaxationRate::setParameters"; +const string WhereGetParameters = "VibrationalRelaxationRate::getParameters"; +const string WhereSetContext = "VibrationalRelaxationRate::setContext"; + +const string ModelConstant = "constant"; +const string ModelMultiState = "multi-state-resolved"; +const string ModelStarikovskiy = "starikovskiy"; +const string ModelCastela = "castela"; + + +// helpers to check the correctness of a chosen vibrational model input data. void requireNoKey(const AnyMap& node, const string& key, const string& model, const string& where) { @@ -176,7 +189,7 @@ void registerVibrationalModelConsistency(const Kinetics& kin, const auto existing = modelByFamily.find(family); if (existing != modelByFamily.end() && existing->second != model) { - throw InputFileError("VibrationalRelaxationRate::setContext", input, + throw InputFileError(WhereSetContext, input, "Inconsistent vibration-model for vibrational family '{}'. " "This family was already registered with model '{}', but the " "current reaction uses model '{}'. A given vibrational family " @@ -201,7 +214,7 @@ void registerVibrationalModelConsistency(const Kinetics& kin, msg << " - " << item.first << ": " << item.second << "\n"; } - warn_user("VibrationalRelaxationRate::setContext", msg.str()); + warn_user(WhereSetContext, msg.str()); s_warnedMixedModels.insert(&kin); } } @@ -211,7 +224,7 @@ string inferRelaxingFamily(const Reaction& rxn) const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); if (vibReactants.empty()) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "A vibrational-relaxation reaction must contain at least one " "vibrational reactant, for example O2(v1), N2(v), or NH3(v1)."); } @@ -226,13 +239,13 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); if (vibReactants.size() != 1) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "vibration-model '{}' expects exactly one vibrational reactant.", model); } if (!vibProducts.empty()) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "vibration-model '{}' describes relaxation to the ground state " "and therefore does not allow vibrational products.", model); } @@ -241,7 +254,7 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, replaceVibrationalSpeciesByGroundState(rxn.reactants); if (!sameComposition(relaxedReactants, rxn.products)) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "Invalid vibrational relaxation stoichiometry for model '{}'. " "Expected a reaction equivalent to X(v) + M => X + M, where the " "collider M is unchanged.", model); @@ -250,7 +263,7 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "vibration-model '{}' expects a bimolecular relaxation reaction " "of the form X(v) + M => X + M.", model); } @@ -258,13 +271,13 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, void validateCastelaReaction(const Reaction& rxn) { - validateSimpleRelaxationToGroundState(rxn, "castela"); + validateSimpleRelaxationToGroundState(rxn, ModelCastela); const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); const string family = vibrationalFamilyName(vibReactants.front()); if (family != "N2(v)") { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "The Castela relaxation model is only valid for N2 vibrational " "relaxation. Found vibrational family '{}'.", family); } @@ -289,7 +302,7 @@ void validateCastelaReaction(const Reaction& rxn) } if (collider != "N2" && collider != "O2" && collider != "O") { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "The Castela relaxation model only supports colliders 'N2', " "'O2', and 'O'. Found collider '{}'.", collider); } @@ -301,7 +314,7 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); if (vibReactants.empty()) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "A detailed vibrational relaxation reaction must contain at least " "one vibrational reactant."); } @@ -313,7 +326,7 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) for (const auto& sp : vibReactants) { const string spFamily = vibrationalFamilyName(sp); if (spFamily != family) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "Invalid detailed vibrational relaxation reaction: all " "vibrational reactants must belong to the same vibrational " "family. Found '{}' and '{}'.", family, spFamily); @@ -323,7 +336,7 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) for (const auto& sp : vibProducts) { const string spFamily = vibrationalFamilyName(sp); if (spFamily != family) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "Invalid detailed vibrational relaxation reaction: all " "vibrational products must belong to the same vibrational " "family. Found '{}' and '{}'.", family, spFamily); @@ -339,7 +352,7 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) replaceVibrationalSpeciesByGroundState(rxn.products); if (!sameComposition(collapsedReactants, collapsedProducts)) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "Invalid detailed vibrational relaxation reaction: replacing all " "vibrational species by their ground-state species does not " "conserve stoichiometry."); @@ -349,7 +362,7 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "Invalid detailed vibrational relaxation reaction: expected a " "bimolecular reaction with two reactant molecules and two product " "molecules."); @@ -471,10 +484,10 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // // The default model is multi-state-resolved. - m_vibration_model = node.getString("vibration-model", "multi-state-resolved"); + m_vibration_model = node.getString("vibration-model", ModelMultiState); if (!node.hasKey("rate-constant")) { - throw InputFileError("VibrationalRelaxationRate::setParameters", node, + throw InputFileError(WhereSetParameters, node, "A vibrational-relaxation reaction requires a 'rate-constant' " "mapping."); } @@ -482,36 +495,36 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, const auto& rate = node["rate-constant"]; if (!rate.is()) { - throw InputFileError("VibrationalRelaxationRate::setParameters", node, + throw InputFileError(WhereSetParameters, node, "The 'rate-constant' field must be a mapping."); } const auto& rateMap = rate.as(); - if (m_vibration_model == "constant") { + if (m_vibration_model == ModelConstant) { // Constant model: // // k(T) = A requireKey(rateMap, m_A_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_b_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, "n", m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_B_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_C_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_D_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_m_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_E_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_z_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_scaling_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); @@ -523,7 +536,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, m_z = 1.0; m_scaling = 1.0; } - else if (m_vibration_model == "multi-state-resolved") { + else if (m_vibration_model == ModelMultiState) { // Detailed VV/VT model: // // k(T) = scaling * A * exp( @@ -534,15 +547,15 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // ) requireKey(rateMap, m_A_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, "n", m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_m_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_E_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_z_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); ArrheniusBase::setParameters(node, rate_units); @@ -554,7 +567,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, m_z = 1.0; m_scaling = rateMap.getDouble(m_scaling_str, 1.0); } - else if (m_vibration_model == "starikovskiy") { + else if (m_vibration_model == ModelStarikovskiy) { // User-facing formula: // // k(T) = A * T^n * exp( @@ -575,11 +588,11 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // m_z = z requireKey(rateMap, m_A_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_b_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_scaling_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); const double n = rateMap.hasKey("n") ? rateMap["n"].asDouble() : 0.0; configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], n); @@ -593,11 +606,11 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, m_scaling = 1.0; if (m_m <= 0.0 || m_z <= 0.0) { - throw InputFileError("VibrationalRelaxationRate::setParameters", node, + throw InputFileError(WhereSetParameters, node, "The Starikovskiy exponents 'm' and 'z' must be positive."); } } - else if (m_vibration_model == "castela") { + else if (m_vibration_model == ModelCastela) { // Castela model: // // Original relaxation time: @@ -621,29 +634,29 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // scaling = 1 requireKey(rateMap, "a", m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireKey(rateMap, "b", m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_A_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, "n", m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, "K", m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_B_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_C_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_D_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_m_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_E_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_z_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); requireNoKey(rateMap, m_scaling_str, m_vibration_model, - "VibrationalRelaxationRate::setParameters"); + WhereSetParameters); m_castela_a = rateMap["a"].asDouble(); m_castela_b = rateMap["b"].asDouble(); @@ -655,7 +668,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, } if (m_referencePressure <= 0.0) { - throw InputFileError("VibrationalRelaxationRate::setParameters", node, + throw InputFileError(WhereSetParameters, node, "Castela reference-pressure must be positive."); } @@ -671,7 +684,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, m_scaling = 1.0; } else { - throw InputFileError("VibrationalRelaxationRate::setParameters", node, + throw InputFileError(WhereSetParameters, node, "Unrecognized vibration-model '{}'. Expected 'multi-state-resolved', " "'starikovskiy', 'castela', or 'constant'.", m_vibration_model); @@ -701,7 +714,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const } }; - if (m_vibration_model == "constant") { + if (m_vibration_model == ModelConstant) { const double tol = 1e-12; if (std::abs(m_b) > tol @@ -710,14 +723,14 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const || std::abs(m_D) > tol || std::abs(m_E) > tol) { - throw InputFileError("VibrationalRelaxationRate::getParameters", node, + throw InputFileError(WhereGetParameters, node, "Cannot serialize this rate as 'constant': the internal " "parameters contain temperature-dependent terms."); } storePreExponentialFactor(rateNode, m_scaling * m_A); } - else if (m_vibration_model == "multi-state-resolved") { + else if (m_vibration_model == ModelMultiState) { storePreExponentialFactor(rateNode, m_A); rateNode[m_b_str] = m_b; @@ -726,7 +739,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const rateNode[m_D_str] = m_D; rateNode[m_scaling_str] = m_scaling; } - else if (m_vibration_model == "starikovskiy") { + else if (m_vibration_model == ModelStarikovskiy) { storePreExponentialFactor(rateNode, m_A); rateNode["n"] = m_b; @@ -737,7 +750,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const rateNode["D"] = m_E; rateNode["z"] = m_z; } - else if (m_vibration_model == "castela") { + else if (m_vibration_model == ModelCastela) { const double tol = 1e-12; if (std::abs(m_b - 1.0) > tol @@ -745,14 +758,14 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const || std::abs(m_E) > tol || std::abs(m_scaling - 1.0) > tol) { - throw InputFileError("VibrationalRelaxationRate::getParameters", node, + throw InputFileError(WhereGetParameters, node, "Cannot serialize this rate as 'castela': the internal " "parameters are not consistent with the Castela form. " "Expected b = 1, D = 0, E = 0, and scaling = 1."); } if (m_referencePressure <= 0.0) { - throw InputFileError("VibrationalRelaxationRate::getParameters", node, + throw InputFileError(WhereGetParameters, node, "Cannot serialize this rate as 'castela': " "reference-pressure must be positive."); } @@ -762,7 +775,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const rateNode[m_reference_pressure_str].setQuantity(m_referencePressure, "Pa"); } else { - throw InputFileError("VibrationalRelaxationRate::getParameters", node, + throw InputFileError(WhereGetParameters, node, "Unrecognized vibration-model '{}'. Expected 'multi-state-resolved', " "'starikovskiy', 'castela', or 'constant'.", m_vibration_model); @@ -787,23 +800,23 @@ double VibrationalRelaxationRate::ddTScaledFromStruct( void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& kin) { if (rxn.reversible) { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "Vibrational relaxation rates do not support reversible " "reactions."); } const string family = inferRelaxingFamily(rxn); - if (m_vibration_model == "constant") { - validateSimpleRelaxationToGroundState(rxn, "constant"); - } else if (m_vibration_model == "castela") { + if (m_vibration_model == ModelConstant) { + validateSimpleRelaxationToGroundState(rxn, ModelConstant); + } else if (m_vibration_model == ModelCastela) { validateCastelaReaction(rxn); - } else if (m_vibration_model == "starikovskiy") { - validateSimpleRelaxationToGroundState(rxn, "starikovskiy"); - } else if (m_vibration_model == "multi-state-resolved") { + } else if (m_vibration_model == ModelStarikovskiy) { + validateSimpleRelaxationToGroundState(rxn, ModelStarikovskiy); + } else if (m_vibration_model == ModelMultiState) { validateDetailedRelaxationReaction(rxn); } else { - throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + throw InputFileError(WhereSetContext, rxn.input, "Unrecognized vibration-model '{}'.", m_vibration_model); } From 624c0fce677ec682915c9ac9ccc7d71e756d02f1 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 22:43:25 +0200 Subject: [PATCH 17/27] make required and forbidden key tests more compact --- src/kinetics/VibrationalRelaxationRate.cpp | 89 +++++++--------------- 1 file changed, 28 insertions(+), 61 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 3243295a325..184dd3603d5 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -53,6 +53,22 @@ void requireKey(const AnyMap& node, const string& key, } } +void requireKeys(const AnyMap& node, const string& model, + const string& where, std::initializer_list keys) +{ + for (const auto& key : keys) { + requireKey(node, key, model, where); + } +} + +void forbidKeys(const AnyMap& node, const string& model, + const string& where, std::initializer_list keys) +{ + for (const auto& key : keys) { + requireNoKey(node, key, model, where); + } +} + // check whether a species is a vibrational excited species // or a vibrational reservoir from its name bool isVibrationalSpecies(const string& name) @@ -505,26 +521,11 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // Constant model: // // k(T) = A - requireKey(rateMap, m_A_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_b_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, "n", m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_B_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_C_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_D_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_m_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_E_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_z_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_scaling_str, m_vibration_model, - WhereSetParameters); + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {m_b_str, "n", m_B_str, m_C_str, m_D_str, m_m_str, + m_E_str, m_z_str, m_scaling_str}); configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); @@ -545,17 +546,8 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // + C * T^(-1/3) // + D * T^(-2/3) // ) - - requireKey(rateMap, m_A_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, "n", m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_m_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_E_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_z_str, m_vibration_model, - WhereSetParameters); + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {"n", m_m_str, m_E_str, m_z_str}); ArrheniusBase::setParameters(node, rate_units); @@ -587,12 +579,8 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // m_E = D // m_z = z - requireKey(rateMap, m_A_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_b_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_scaling_str, m_vibration_model, - WhereSetParameters); + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {m_b_str, m_scaling_str}); const double n = rateMap.hasKey("n") ? rateMap["n"].asDouble() : 0.0; configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], n); @@ -633,30 +621,9 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // E = 0 // scaling = 1 - requireKey(rateMap, "a", m_vibration_model, - WhereSetParameters); - requireKey(rateMap, "b", m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_A_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, "n", m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, "K", m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_B_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_C_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_D_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_m_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_E_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_z_str, m_vibration_model, - WhereSetParameters); - requireNoKey(rateMap, m_scaling_str, m_vibration_model, - WhereSetParameters); + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {"a", "b"}); + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str, "n", "K", + m_B_str, m_C_str, m_D_str, m_m_str, m_E_str, m_z_str, m_scaling_str}); m_castela_a = rateMap["a"].asDouble(); m_castela_b = rateMap["b"].asDouble(); From 32cb13bbf08f9a6a0963d62a4d25d0dfe5aff08f Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 23:04:40 +0200 Subject: [PATCH 18/27] make data input in setParameters more readable and compact --- .../kinetics/VibrationalRelaxationRate.h | 4 ++ src/kinetics/VibrationalRelaxationRate.cpp | 59 +++++++------------ 2 files changed, 24 insertions(+), 39 deletions(-) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index 825dad50a6a..45e0aa12baf 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -271,6 +271,10 @@ class VibrationalRelaxationRate : public ArrheniusBase */ void configureBaseFromYamlA(const AnyMap& node, const UnitStack& rate_units, const AnyValue& A, double b); + + //! Sets parameters + void setGenericParameters(double B, double C, double D, double m, + double E, double z, double scaling); }; } diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 184dd3603d5..9032b9cb876 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -528,14 +528,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, m_E_str, m_z_str, m_scaling_str}); configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); - - m_B = 0.0; - m_C = 0.0; - m_D = 0.0; - m_m = 2.0 / 3.0; - m_E = 0.0; - m_z = 1.0; - m_scaling = 1.0; + setGenericParameters(0.0, 0.0, 0.0, 2.3/3.0, 0.0, 1.0, 1.0); } else if (m_vibration_model == ModelMultiState) { // Detailed VV/VT model: @@ -550,14 +543,8 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {"n", m_m_str, m_E_str, m_z_str}); ArrheniusBase::setParameters(node, rate_units); - - m_B = rateMap.getDouble(m_B_str, 0.0); - m_C = rateMap.getDouble(m_C_str, 0.0); - m_D = rateMap.getDouble(m_D_str, 0.0); - m_m = 2.0 / 3.0; - m_E = 0.0; - m_z = 1.0; - m_scaling = rateMap.getDouble(m_scaling_str, 1.0); + setGenericParameters(rateMap.getDouble(m_B_str, 0.0), rateMap.getDouble(m_C_str, 0.0), + rateMap.getDouble(m_D_str, 0.0), 2.0 / 3.0, 0.0, 1.0, rateMap.getDouble(m_scaling_str, 1.0)); } else if (m_vibration_model == ModelStarikovskiy) { // User-facing formula: @@ -582,16 +569,9 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {m_b_str, m_scaling_str}); - const double n = rateMap.hasKey("n") ? rateMap["n"].asDouble() : 0.0; - configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], n); - - m_B = rateMap.getDouble("K", 0.0); - m_C = rateMap.getDouble("B", 0.0); - m_D = rateMap.getDouble("C", 0.0); - m_m = rateMap.getDouble("m", 1.0); - m_E = rateMap.getDouble("D", 0.0); - m_z = rateMap.getDouble("z", 1.0); - m_scaling = 1.0; + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], rateMap.getDouble("n", 0.0)); + setGenericParameters(rateMap.getDouble("K", 0.0), rateMap.getDouble("B", 0.0), rateMap.getDouble("C", 0.0), + rateMap.getDouble("m", 1.0), rateMap.getDouble("D", 0.0), rateMap.getDouble("z", 1.0), 1.0); if (m_m <= 0.0 || m_z <= 0.0) { throw InputFileError(WhereSetParameters, node, @@ -625,9 +605,6 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str, "n", "K", m_B_str, m_C_str, m_D_str, m_m_str, m_E_str, m_z_str, m_scaling_str}); - m_castela_a = rateMap["a"].asDouble(); - m_castela_b = rateMap["b"].asDouble(); - if (rateMap.hasKey(m_reference_pressure_str)) { m_referencePressure = rateMap.convert(m_reference_pressure_str, "Pa"); } else { @@ -639,16 +616,8 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, "Castela reference-pressure must be positive."); } - configureBaseFromInternalA( - node, rate_units, GasConstant / m_referencePressure, 1.0); - - m_B = 18.42 + m_castela_a * m_castela_b; - m_C = -m_castela_a; - m_D = 0.0; - m_m = 2.0 / 3.0; - m_E = 0.0; - m_z = 1.0; - m_scaling = 1.0; + configureBaseFromInternalA(node, rate_units, GasConstant / m_referencePressure, 1.0); + setGenericParameters(18.42 + rateMap["a"].asDouble() * rateMap["b"].asDouble(), -rateMap["a"].asDouble(), 0.0, 2.0 / 3.0, 0.0, 1.0, 1.0); } else { throw InputFileError(WhereSetParameters, node, @@ -791,4 +760,16 @@ void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& kin, family, m_vibration_model, rxn.input); } +void VibrationalRelaxationRate::setGenericParameters( + double B, double C, double D, double m, double E, double z, double scaling) +{ + m_B = B; + m_C = C; + m_D = D; + m_m = m; + m_E = E; + m_z = z; + m_scaling = scaling; +} + } // namespace Cantera \ No newline at end of file From 133d934aea635eec4a7860082badb239f0dbbbf3 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 23:21:39 +0200 Subject: [PATCH 19/27] make sameComposition more compact and readable --- src/kinetics/VibrationalRelaxationRate.cpp | 30 +++++----------------- 1 file changed, 6 insertions(+), 24 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 9032b9cb876..83e8dc7aa7b 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -127,33 +127,15 @@ Composition replaceVibrationalSpeciesByGroundState(const Composition& comp) return out; } -bool sameComposition(const Composition& a, const Composition& b, - double tol = 1e-12) +bool sameComposition(Composition diff, const Composition& b, + double tol = Tiny) { - std::set names; - - for (const auto& item : a) { - names.insert(item.first); - } - for (const auto& item : b) { - names.insert(item.first); + for (const auto& [name, value] : b) { + diff[name] -= value; } - for (const auto& name : names) { - double av = 0.0; - double bv = 0.0; - - const auto ait = a.find(name); - if (ait != a.end()) { - av = ait->second; - } - - const auto bit = b.find(name); - if (bit != b.end()) { - bv = bit->second; - } - - if (std::abs(av - bv) > tol) { + for (const auto& [name, value] : diff) { + if (std::abs(value) > tol) { return false; } } From 06a051b62f0deb0e65c7e5b6c929c257180981d9 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 23:36:39 +0200 Subject: [PATCH 20/27] correct typos and finish corrections that were forgotten during previous commits --- src/kinetics/VibrationalRelaxationRate.cpp | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 83e8dc7aa7b..9ef68b6dfff 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -258,8 +258,8 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, "collider M is unchanged.", model); } - if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 - || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) + if (std::abs(compositionSum(rxn.reactants) - 2.0) > Tiny + || std::abs(compositionSum(rxn.products) - 2.0) > Tiny) { throw InputFileError(WhereSetContext, rxn.input, "vibration-model '{}' expects a bimolecular relaxation reaction " @@ -291,10 +291,10 @@ void validateCastelaReaction(const Reaction& rxn) const double value = item.second; if (name == "N2") { - if (std::abs(value - 2.0) < 1e-12) { + if (std::abs(value - 2.0) < Tiny) { collider = "N2"; } - } else if (std::abs(value - 1.0) < 1e-12) { + } else if (std::abs(value - 1.0) < Tiny) { collider = name; } } @@ -357,8 +357,8 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) } // The current detailed VV/VT formulation is bimolecular. - if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 - || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) + if (std::abs(compositionSum(rxn.reactants) - 2.0) > Tiny + || std::abs(compositionSum(rxn.products) - 2.0) > Tiny) { throw InputFileError(WhereSetContext, rxn.input, "Invalid detailed vibrational relaxation reaction: expected a " @@ -510,7 +510,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, m_E_str, m_z_str, m_scaling_str}); configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); - setGenericParameters(0.0, 0.0, 0.0, 2.3/3.0, 0.0, 1.0, 1.0); + setGenericParameters(0.0, 0.0, 0.0, 2.0/3.0, 0.0, 1.0, 1.0); } else if (m_vibration_model == ModelMultiState) { // Detailed VV/VT model: @@ -599,7 +599,9 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, } configureBaseFromInternalA(node, rate_units, GasConstant / m_referencePressure, 1.0); - setGenericParameters(18.42 + rateMap["a"].asDouble() * rateMap["b"].asDouble(), -rateMap["a"].asDouble(), 0.0, 2.0 / 3.0, 0.0, 1.0, 1.0); + m_castela_a = rateMap["a"].asDouble(); + m_castela_b = rateMap["b"].asDouble(); + setGenericParameters(18.42 + m_castela_a * m_castela_b, -m_castela_a, 0.0, 2.0 / 3.0, 0.0, 1.0, 1.0); } else { throw InputFileError(WhereSetParameters, node, @@ -633,7 +635,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const }; if (m_vibration_model == ModelConstant) { - const double tol = 1e-12; + const double tol = Tiny; if (std::abs(m_b) > tol || std::abs(m_B) > tol @@ -669,7 +671,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const rateNode["z"] = m_z; } else if (m_vibration_model == ModelCastela) { - const double tol = 1e-12; + const double tol = Tiny; if (std::abs(m_b - 1.0) > tol || std::abs(m_D) > tol From 72d98c6782d6e992503892604a8ac4076cbd78b4 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 23:45:18 +0200 Subject: [PATCH 21/27] Factor vibrational relaxation rate-constant map parsing and correct Tiny declaration error --- src/kinetics/VibrationalRelaxationRate.cpp | 36 ++++++++++++---------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 9ef68b6dfff..514292b11d7 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -33,6 +33,23 @@ const string ModelMultiState = "multi-state-resolved"; const string ModelStarikovskiy = "starikovskiy"; const string ModelCastela = "castela"; +const AnyMap& getRateConstantMap(const AnyMap& node) +{ + if (!node.hasKey("rate-constant")) { + throw InputFileError(WhereSetParameters, node, + "A vibrational-relaxation reaction requires a 'rate-constant' " + "mapping."); + } + + const auto& rate = node["rate-constant"]; + + if (!rate.is()) { + throw InputFileError(WhereSetParameters, node, + "The 'rate-constant' field must be a mapping."); + } + + return rate.as(); +} // helpers to check the correctness of a chosen vibrational model input data. void requireNoKey(const AnyMap& node, const string& key, @@ -484,20 +501,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, m_vibration_model = node.getString("vibration-model", ModelMultiState); - if (!node.hasKey("rate-constant")) { - throw InputFileError(WhereSetParameters, node, - "A vibrational-relaxation reaction requires a 'rate-constant' " - "mapping."); - } - - const auto& rate = node["rate-constant"]; - - if (!rate.is()) { - throw InputFileError(WhereSetParameters, node, - "The 'rate-constant' field must be a mapping."); - } - - const auto& rateMap = rate.as(); + const auto& rateMap = getRateConstantMap(node); if (m_vibration_model == ModelConstant) { // Constant model: @@ -635,7 +639,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const }; if (m_vibration_model == ModelConstant) { - const double tol = Tiny; + const double tol = 1e-12; if (std::abs(m_b) > tol || std::abs(m_B) > tol @@ -671,7 +675,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const rateNode["z"] = m_z; } else if (m_vibration_model == ModelCastela) { - const double tol = Tiny; + const double tol = 1e-12; if (std::abs(m_b - 1.0) > tol || std::abs(m_D) > tol From 5b2e7353819ebdc37caf658cc26e14fe12b54b84 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 9 Jul 2026 23:55:46 +0200 Subject: [PATCH 22/27] Split vibrational relaxation parameter treatment by model and use switch / case --- .../kinetics/VibrationalRelaxationRate.h | 16 + src/kinetics/VibrationalRelaxationRate.cpp | 294 +++++++++++------- 2 files changed, 196 insertions(+), 114 deletions(-) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index 45e0aa12baf..f3be755c7c1 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -275,6 +275,22 @@ class VibrationalRelaxationRate : public ArrheniusBase //! Sets parameters void setGenericParameters(double B, double C, double D, double m, double E, double z, double scaling); + + //! Sub-function of setParameters relative to the 'constant' model + void setConstantParameters(const AnyMap& node, const AnyMap& rateMap, + const UnitStack& rate_units); + + //! Sub-function of setParameters relative to the 'multi-state-resolved' model + void setMultiStateParameters(const AnyMap& node, const AnyMap& rateMap, + const UnitStack& rate_units); + + //! Sub-function of setParameters relative to the 'starikovskiy' model + void setStarikovskiyParameters(const AnyMap& node, const AnyMap& rateMap, + const UnitStack& rate_units); + + //! Sub-function of setParameters relative to the 'castela' model + void setCastelaParameters(const AnyMap& node, const AnyMap& rateMap, + const UnitStack& rate_units); }; } diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 514292b11d7..d0d241e5033 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -33,6 +33,31 @@ const string ModelMultiState = "multi-state-resolved"; const string ModelStarikovskiy = "starikovskiy"; const string ModelCastela = "castela"; +enum class VibModel { + Constant, + MultiStateResolved, + Starikovskiy, + Castela +}; + +VibModel parseVibrationModel(const string& model, const AnyMap& input, + const string& where) +{ + if (model == ModelConstant) { + return VibModel::Constant; + } else if (model == ModelMultiState) { + return VibModel::MultiStateResolved; + } else if (model == ModelStarikovskiy) { + return VibModel::Starikovskiy; + } else if (model == ModelCastela) { + return VibModel::Castela; + } + + throw InputFileError(where, input, + "Unrecognized vibration-model '{}'. Expected 'multi-state-resolved', " + "'starikovskiy', 'castela', or 'constant'.", model); +} + const AnyMap& getRateConstantMap(const AnyMap& node) { if (!node.hasKey("rate-constant")) { @@ -493,126 +518,165 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // vibration-model: starikovskiy // vibration-model: castela // - // This is done in a spriti of class adaptability: - // Any user wishing to develop more advanced vibration models - // is welcome to add them in this class. - // // The default model is multi-state-resolved. m_vibration_model = node.getString("vibration-model", ModelMultiState); - const auto& rateMap = getRateConstantMap(node); - if (m_vibration_model == ModelConstant) { - // Constant model: - // - // k(T) = A - requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + switch (parseVibrationModel(m_vibration_model, node, WhereSetParameters)) { + case VibModel::Constant: + setConstantParameters(node, rateMap, rate_units); + break; + case VibModel::MultiStateResolved: + setMultiStateParameters(node, rateMap, rate_units); + break; + case VibModel::Starikovskiy: + setStarikovskiyParameters(node, rateMap, rate_units); + break; + case VibModel::Castela: + setCastelaParameters(node, rateMap, rate_units); + break; + } +} - forbidKeys(rateMap, m_vibration_model, WhereSetParameters, - {m_b_str, "n", m_B_str, m_C_str, m_D_str, m_m_str, - m_E_str, m_z_str, m_scaling_str}); +void VibrationalRelaxationRate::setConstantParameters( + const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units) +{ + // Constant model: + // + // k(T) = A + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); - configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); - setGenericParameters(0.0, 0.0, 0.0, 2.0/3.0, 0.0, 1.0, 1.0); - } - else if (m_vibration_model == ModelMultiState) { - // Detailed VV/VT model: - // - // k(T) = scaling * A * exp( - // b * log(T) - // + B - // + C * T^(-1/3) - // + D * T^(-2/3) - // ) - requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); - forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {"n", m_m_str, m_E_str, m_z_str}); - - ArrheniusBase::setParameters(node, rate_units); - setGenericParameters(rateMap.getDouble(m_B_str, 0.0), rateMap.getDouble(m_C_str, 0.0), - rateMap.getDouble(m_D_str, 0.0), 2.0 / 3.0, 0.0, 1.0, rateMap.getDouble(m_scaling_str, 1.0)); - } - else if (m_vibration_model == ModelStarikovskiy) { - // User-facing formula: - // - // k(T) = A * T^n * exp( - // K - // + B * T^(-1/3) - // + C * T^(-m) - // + D * T^(-z) - // ) - // - // Internal mapping: - // - // m_b = n - // m_B = K - // m_C = B - // m_D = C - // m_m = m - // m_E = D - // m_z = z - - requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); - forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {m_b_str, m_scaling_str}); - - configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], rateMap.getDouble("n", 0.0)); - setGenericParameters(rateMap.getDouble("K", 0.0), rateMap.getDouble("B", 0.0), rateMap.getDouble("C", 0.0), - rateMap.getDouble("m", 1.0), rateMap.getDouble("D", 0.0), rateMap.getDouble("z", 1.0), 1.0); - - if (m_m <= 0.0 || m_z <= 0.0) { - throw InputFileError(WhereSetParameters, node, - "The Starikovskiy exponents 'm' and 'z' must be positive."); - } + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {m_b_str, "n", m_B_str, m_C_str, m_D_str, m_m_str, + m_E_str, m_z_str, m_scaling_str}); + + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); + setGenericParameters(0.0, 0.0, 0.0, 2.0 / 3.0, 0.0, 1.0, 1.0); +} + +void VibrationalRelaxationRate::setMultiStateParameters( + const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units) +{ + // Detailed VV/VT model: + // + // k(T) = scaling * A * exp( + // b * log(T) + // + B + // + C * T^(-1/3) + // + D * T^(-2/3) + // ) + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {"n", m_m_str, m_E_str, m_z_str}); + + ArrheniusBase::setParameters(node, rate_units); + + setGenericParameters( + rateMap.getDouble(m_B_str, 0.0), + rateMap.getDouble(m_C_str, 0.0), + rateMap.getDouble(m_D_str, 0.0), + 2.0 / 3.0, + 0.0, + 1.0, + rateMap.getDouble(m_scaling_str, 1.0)); +} + +void VibrationalRelaxationRate::setStarikovskiyParameters( + const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units) +{ + // User-facing formula: + // + // k(T) = A * T^n * exp( + // K + // + B * T^(-1/3) + // + C * T^(-m) + // + D * T^(-z) + // ) + // + // B, C, and D are signed coefficients read directly from YAML. + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {m_b_str, m_scaling_str}); + + const double m = rateMap.getDouble("m", 1.0); + const double z = rateMap.getDouble("z", 1.0); + + if (m <= 0.0 || z <= 0.0) { + throw InputFileError(WhereSetParameters, node, + "The Starikovskiy exponents 'm' and 'z' must be positive."); } - else if (m_vibration_model == ModelCastela) { - // Castela model: - // - // Original relaxation time: - // - // tau_k = p0 / p_k - // * exp[a_k * (T^(-1/3) - b_k) - 18.42] - // - // Equivalent bimolecular rate coefficient: - // - // k_k(T) = R T / p0 - // * exp[18.42 + a_k b_k - a_k T^(-1/3)] - // - // Internal mapping: - // - // A = R / p0 - // b = 1 - // B = 18.42 + a_k b_k - // C = -a_k - // D = 0 - // E = 0 - // scaling = 1 - - requireKeys(rateMap, m_vibration_model, WhereSetParameters, {"a", "b"}); - forbidKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str, "n", "K", - m_B_str, m_C_str, m_D_str, m_m_str, m_E_str, m_z_str, m_scaling_str}); - - if (rateMap.hasKey(m_reference_pressure_str)) { - m_referencePressure = rateMap.convert(m_reference_pressure_str, "Pa"); - } else { - m_referencePressure = OneAtm; - } - if (m_referencePressure <= 0.0) { - throw InputFileError(WhereSetParameters, node, - "Castela reference-pressure must be positive."); - } + configureBaseFromYamlA( + node, rate_units, rateMap[m_A_str], rateMap.getDouble("n", 0.0)); - configureBaseFromInternalA(node, rate_units, GasConstant / m_referencePressure, 1.0); - m_castela_a = rateMap["a"].asDouble(); - m_castela_b = rateMap["b"].asDouble(); - setGenericParameters(18.42 + m_castela_a * m_castela_b, -m_castela_a, 0.0, 2.0 / 3.0, 0.0, 1.0, 1.0); + setGenericParameters( + rateMap.getDouble("K", 0.0), + rateMap.getDouble("B", 0.0), + rateMap.getDouble("C", 0.0), + m, + rateMap.getDouble("D", 0.0), + z, + 1.0); +} + +void VibrationalRelaxationRate::setCastelaParameters( + const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units) +{ + // Castela model: + // + // Original relaxation time: + // + // tau_k = p0 / p_k + // * exp[a_k * (T^(-1/3) - b_k) - 18.42] + // + // Equivalent bimolecular rate coefficient: + // + // k_k(T) = R T / p0 + // * exp[18.42 + a_k b_k - a_k T^(-1/3)] + // + // Internal mapping: + // + // A = R / p0 + // b = 1 + // B = 18.42 + a_k b_k + // C = -a_k + // D = 0 + // E = 0 + // scaling = 1 + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {"a", "b"}); + + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {m_A_str, "n", "K", m_B_str, m_C_str, m_D_str, m_m_str, + m_E_str, m_z_str, m_scaling_str}); + + m_castela_a = rateMap["a"].asDouble(); + m_castela_b = rateMap["b"].asDouble(); + + if (rateMap.hasKey(m_reference_pressure_str)) { + m_referencePressure = rateMap.convert(m_reference_pressure_str, "Pa"); + } else { + m_referencePressure = OneAtm; } - else { + + if (m_referencePressure <= 0.0) { throw InputFileError(WhereSetParameters, node, - "Unrecognized vibration-model '{}'. Expected 'multi-state-resolved', " - "'starikovskiy', 'castela', or 'constant'.", - m_vibration_model); + "Castela reference-pressure must be positive."); } + + configureBaseFromInternalA( + node, rate_units, GasConstant / m_referencePressure, 1.0); + + setGenericParameters( + 18.42 + m_castela_a * m_castela_b, + -m_castela_a, + 0.0, + 2.0 / 3.0, + 0.0, + 1.0, + 1.0); } void VibrationalRelaxationRate::getParameters(AnyMap& node) const @@ -731,17 +795,19 @@ void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& const string family = inferRelaxingFamily(rxn); - if (m_vibration_model == ModelConstant) { + switch (parseVibrationModel(m_vibration_model, rxn.input, WhereSetContext)) { + case VibModel::Constant: validateSimpleRelaxationToGroundState(rxn, ModelConstant); - } else if (m_vibration_model == ModelCastela) { + break; + case VibModel::Castela: validateCastelaReaction(rxn); - } else if (m_vibration_model == ModelStarikovskiy) { + break; + case VibModel::Starikovskiy: validateSimpleRelaxationToGroundState(rxn, ModelStarikovskiy); - } else if (m_vibration_model == ModelMultiState) { + break; + case VibModel::MultiStateResolved: validateDetailedRelaxationReaction(rxn); - } else { - throw InputFileError(WhereSetContext, rxn.input, - "Unrecognized vibration-model '{}'.", m_vibration_model); + break; } registerVibrationalModelConsistency( From e773816a67dc0896810a3e0daa9ae73a1d21ea47 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 10 Jul 2026 00:12:14 +0200 Subject: [PATCH 23/27] split getParameter by model as well --- .../kinetics/VibrationalRelaxationRate.h | 14 ++ src/kinetics/VibrationalRelaxationRate.cpp | 176 ++++++++++-------- 2 files changed, 111 insertions(+), 79 deletions(-) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index f3be755c7c1..5065ac72eda 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -291,6 +291,20 @@ class VibrationalRelaxationRate : public ArrheniusBase //! Sub-function of setParameters relative to the 'castela' model void setCastelaParameters(const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units); + + void storePreExponentialFactor(AnyMap& target, double A) const; + + //! Sub-function of getParameters relative to the 'constant' model + void getConstantParameters(AnyMap& node, AnyMap& rateNode) const; + + //! Sub-function of getParameters relative to the 'multi-state-resolved' model + void getMultiStateParameters(AnyMap& rateNode) const; + + //! Sub-function of getParameters relative to the 'starikovskiy' model + void getStarikovskiyParameters(AnyMap& rateNode) const; + + //! Sub-function of getParameters relative to the 'castela' model + void getCastelaParameters(AnyMap& node, AnyMap& rateNode) const; }; } diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index d0d241e5033..1c7e3d9eb51 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -22,7 +22,7 @@ namespace Cantera namespace { -constexpr double Tiny = 1e-12; +constexpr double VibTolerance = 1e-12; const string WhereSetParameters = "VibrationalRelaxationRate::setParameters"; const string WhereGetParameters = "VibrationalRelaxationRate::getParameters"; @@ -170,7 +170,7 @@ Composition replaceVibrationalSpeciesByGroundState(const Composition& comp) } bool sameComposition(Composition diff, const Composition& b, - double tol = Tiny) + double tol = VibTolerance) { for (const auto& [name, value] : b) { diff[name] -= value; @@ -300,8 +300,8 @@ void validateSimpleRelaxationToGroundState(const Reaction& rxn, "collider M is unchanged.", model); } - if (std::abs(compositionSum(rxn.reactants) - 2.0) > Tiny - || std::abs(compositionSum(rxn.products) - 2.0) > Tiny) + if (std::abs(compositionSum(rxn.reactants) - 2.0) > VibTolerance + || std::abs(compositionSum(rxn.products) - 2.0) > VibTolerance) { throw InputFileError(WhereSetContext, rxn.input, "vibration-model '{}' expects a bimolecular relaxation reaction " @@ -333,10 +333,10 @@ void validateCastelaReaction(const Reaction& rxn) const double value = item.second; if (name == "N2") { - if (std::abs(value - 2.0) < Tiny) { + if (std::abs(value - 2.0) < VibTolerance) { collider = "N2"; } - } else if (std::abs(value - 1.0) < Tiny) { + } else if (std::abs(value - 1.0) < VibTolerance) { collider = name; } } @@ -399,8 +399,8 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) } // The current detailed VV/VT formulation is bimolecular. - if (std::abs(compositionSum(rxn.reactants) - 2.0) > Tiny - || std::abs(compositionSum(rxn.products) - 2.0) > Tiny) + if (std::abs(compositionSum(rxn.reactants) - 2.0) > VibTolerance + || std::abs(compositionSum(rxn.products) - 2.0) > VibTolerance) { throw InputFileError(WhereSetContext, rxn.input, "Invalid detailed vibrational relaxation reaction: expected a " @@ -693,84 +693,102 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const AnyMap rateNode; - auto storePreExponentialFactor = [&](AnyMap& target, double A) { - if (conversionUnits().factor() != 0.0) { - target[m_A_str].setQuantity(A, conversionUnits()); - } else { - target[m_A_str] = A; - target["__unconvertible__"] = true; - } - }; - - if (m_vibration_model == ModelConstant) { - const double tol = 1e-12; - - if (std::abs(m_b) > tol - || std::abs(m_B) > tol - || std::abs(m_C) > tol - || std::abs(m_D) > tol - || std::abs(m_E) > tol) - { - throw InputFileError(WhereGetParameters, node, - "Cannot serialize this rate as 'constant': the internal " - "parameters contain temperature-dependent terms."); - } + switch (parseVibrationModel(m_vibration_model, node, WhereGetParameters)) { + case VibModel::Constant: + getConstantParameters(node, rateNode); + break; + case VibModel::MultiStateResolved: + getMultiStateParameters(rateNode); + break; + case VibModel::Starikovskiy: + getStarikovskiyParameters(rateNode); + break; + case VibModel::Castela: + getCastelaParameters(node, rateNode); + break; + } - storePreExponentialFactor(rateNode, m_scaling * m_A); - } - else if (m_vibration_model == ModelMultiState) { - storePreExponentialFactor(rateNode, m_A); - - rateNode[m_b_str] = m_b; - rateNode[m_B_str] = m_B; - rateNode[m_C_str] = m_C; - rateNode[m_D_str] = m_D; - rateNode[m_scaling_str] = m_scaling; - } - else if (m_vibration_model == ModelStarikovskiy) { - storePreExponentialFactor(rateNode, m_A); - - rateNode["n"] = m_b; - rateNode["K"] = m_B; - rateNode["B"] = m_C; - rateNode["C"] = m_D; - rateNode["m"] = m_m; - rateNode["D"] = m_E; - rateNode["z"] = m_z; - } - else if (m_vibration_model == ModelCastela) { - const double tol = 1e-12; - - if (std::abs(m_b - 1.0) > tol - || std::abs(m_D) > tol - || std::abs(m_E) > tol - || std::abs(m_scaling - 1.0) > tol) - { - throw InputFileError(WhereGetParameters, node, - "Cannot serialize this rate as 'castela': the internal " - "parameters are not consistent with the Castela form. " - "Expected b = 1, D = 0, E = 0, and scaling = 1."); - } + rateNode.setFlowStyle(); + node["rate-constant"] = std::move(rateNode); +} - if (m_referencePressure <= 0.0) { - throw InputFileError(WhereGetParameters, node, - "Cannot serialize this rate as 'castela': " - "reference-pressure must be positive."); - } +void VibrationalRelaxationRate::storePreExponentialFactor( + AnyMap& target, double A) const +{ + if (conversionUnits().factor() != 0.0) { + target[m_A_str].setQuantity(A, conversionUnits()); + } else { + target[m_A_str] = A; + target["__unconvertible__"] = true; + } +} - rateNode["a"] = m_castela_a; - rateNode["b"] = m_castela_b; - rateNode[m_reference_pressure_str].setQuantity(m_referencePressure, "Pa"); +void VibrationalRelaxationRate::getConstantParameters( + AnyMap& node, AnyMap& rateNode) const +{ + if (std::abs(m_b) > VibTolerance + || std::abs(m_B) > VibTolerance + || std::abs(m_C) > VibTolerance + || std::abs(m_D) > VibTolerance + || std::abs(m_E) > VibTolerance) + { + throw InputFileError(WhereGetParameters, node, + "Cannot serialize this rate as 'constant': the internal " + "parameters contain temperature-dependent terms."); } - else { + + storePreExponentialFactor(rateNode, m_scaling * m_A); +} + +void VibrationalRelaxationRate::getMultiStateParameters( + AnyMap& rateNode) const +{ + storePreExponentialFactor(rateNode, m_A); + + rateNode[m_b_str] = m_b; + rateNode[m_B_str] = m_B; + rateNode[m_C_str] = m_C; + rateNode[m_D_str] = m_D; + rateNode[m_scaling_str] = m_scaling; +} + +void VibrationalRelaxationRate::getStarikovskiyParameters( + AnyMap& rateNode) const +{ + storePreExponentialFactor(rateNode, m_A); + + rateNode["n"] = m_b; + rateNode["K"] = m_B; + rateNode["B"] = m_C; + rateNode["C"] = m_D; + rateNode["m"] = m_m; + rateNode["D"] = m_E; + rateNode["z"] = m_z; +} + +void VibrationalRelaxationRate::getCastelaParameters( + AnyMap& node, AnyMap& rateNode) const +{ + if (std::abs(m_b - 1.0) > VibTolerance + || std::abs(m_D) > VibTolerance + || std::abs(m_E) > VibTolerance + || std::abs(m_scaling - 1.0) > VibTolerance) + { throw InputFileError(WhereGetParameters, node, - "Unrecognized vibration-model '{}'. Expected 'multi-state-resolved', " - "'starikovskiy', 'castela', or 'constant'.", - m_vibration_model); + "Cannot serialize this rate as 'castela': the internal " + "parameters are not consistent with the Castela form. " + "Expected b = 1, D = 0, E = 0, and scaling = 1."); } - rateNode.setFlowStyle(); - node["rate-constant"] = std::move(rateNode); + if (m_referencePressure <= 0.0) { + throw InputFileError(WhereGetParameters, node, + "Cannot serialize this rate as 'castela': " + "reference-pressure must be positive."); + } + + rateNode["a"] = m_castela_a; + rateNode["b"] = m_castela_b; + rateNode[m_reference_pressure_str].setQuantity(m_referencePressure, "Pa"); } double VibrationalRelaxationRate::ddTScaledFromStruct( From 162e01e09ad7c0b9312a7ca63ed251c2bcb747c5 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 10 Jul 2026 00:15:37 +0200 Subject: [PATCH 24/27] validateDetailedRelaxationReaction made a little more compact and readable --- src/kinetics/VibrationalRelaxationRate.cpp | 37 +++++++++++----------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 1c7e3d9eb51..2016b262772 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -348,6 +348,22 @@ void validateCastelaReaction(const Reaction& rxn) } } +void validateSameVibrationalFamily(const std::vector& species, + const string& family, + const AnyMap& input, + const string& role) +{ + for (const auto& sp : species) { + const string spFamily = vibrationalFamilyName(sp); + if (spFamily != family) { + throw InputFileError(WhereSetContext, input, + "Invalid detailed vibrational relaxation reaction: all " + "vibrational {} must belong to the same vibrational " + "family. Found '{}' and '{}'.", role, family, spFamily); + } + } +} + void validateDetailedRelaxationReaction(const Reaction& rxn) { const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); @@ -363,25 +379,8 @@ void validateDetailedRelaxationReaction(const Reaction& rxn) // to belong to the same "vibrational family". const string family = vibrationalFamilyName(vibReactants.front()); - for (const auto& sp : vibReactants) { - const string spFamily = vibrationalFamilyName(sp); - if (spFamily != family) { - throw InputFileError(WhereSetContext, rxn.input, - "Invalid detailed vibrational relaxation reaction: all " - "vibrational reactants must belong to the same vibrational " - "family. Found '{}' and '{}'.", family, spFamily); - } - } - - for (const auto& sp : vibProducts) { - const string spFamily = vibrationalFamilyName(sp); - if (spFamily != family) { - throw InputFileError(WhereSetContext, rxn.input, - "Invalid detailed vibrational relaxation reaction: all " - "vibrational products must belong to the same vibrational " - "family. Found '{}' and '{}'.", family, spFamily); - } - } + validateSameVibrationalFamily(vibReactants, family, rxn.input, "reactants"); + validateSameVibrationalFamily(vibProducts, family, rxn.input, "products"); // The reaction must conserve the ground-state composition once all // vibrational labels are removed. From cc10c0c248083b885e1560a9054b6941acfe9daa Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 10 Jul 2026 00:50:30 +0200 Subject: [PATCH 25/27] update reactions documentation for the vibration-relaxation reaction type --- doc/sphinx/yaml/reactions.md | 280 +++++++++++++++++++++++++++++++++++ 1 file changed, 280 insertions(+) diff --git a/doc/sphinx/yaml/reactions.md b/doc/sphinx/yaml/reactions.md index 44496a4a5ae..84b51686572 100644 --- a/doc/sphinx/yaml/reactions.md +++ b/doc/sphinx/yaml/reactions.md @@ -18,6 +18,7 @@ The fields common to all `reaction` entries are: - [`Blowers-Masel`](sec-yaml-Blowers-Masel) - [`two-temperature-plasma`](sec-yaml-two-temperature-plasma) - [`electron-collision-plasma`](sec-yaml-electron-collision-plasma) + - [`vibrational-relaxation`](sec-yaml-vibrational-relaxation) - [`electron-collisions`](sec-yaml-electron-collisions) - [`falloff`](sec-yaml-falloff) - [`chemically-activated`](sec-yaml-chemically-activated) @@ -135,6 +136,240 @@ or a corresponding four-element list. The following are equivalent: {A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol, b-gas: 0} [17283, -3.1, -5820 J/mol, 1081 J/mol] ``` +(sec-vibrational-relaxation-constant-rate)= + +### Vibrational relaxation: constant + +The `constant` vibrational relaxation rate expression is + +$$ +k_f = A +$$ + +where $A$ is a temperature-independent rate coefficient. + +The rate parameters are specified as a mapping with fields: + +`A` +: The pre-exponential factor. The units are the standard Cantera rate coefficient +units for the reaction order and are converted using the unit system. + +No other rate parameters are allowed for this model. In particular, `b`, `n`, `B`, +`C`, `D`, `m`, `E`, `z`, and `scaling` are not permitted. + +This model describes relaxation of one vibrationally excited reactant to its ground +state by collision with an unchanged collider: + +```yaml +- equation: N2(v) + N2 => N2 + N2 + type: vibrational-relaxation + vibration-model: constant + rate-constant: {A: 1.0e2 cm^3/mol/s} +``` + +(sec-vibrational-relaxation-multi-state-resolved-rate)= + +### Vibrational relaxation: multi-state-resolved + +The `multi-state-resolved` vibrational relaxation rate expression is + +$$ +k_f = +\mathrm{scaling} , A +\exp \left( +b \ln T ++ B ++ C T^{-1/3} ++ D T^{-2/3} +\right) +$$ + +where $T$ is the gas temperature in K. + +The rate parameters are specified as a mapping with fields: + +`A` +: The pre-exponential factor. The units are the standard Cantera rate coefficient +units for the reaction order and are converted using the unit system. + +`b` +: The dimensionless temperature exponent. Defaults to 0.0. + +`B` +: Dimensionless constant term in the exponential. Defaults to 0.0. + +`C` +: Coefficient multiplying $T^{-1/3}$. This value is interpreted assuming that +$T$ is in K. Defaults to 0.0. + +`D` +: Coefficient multiplying $T^{-2/3}$. This value is interpreted assuming that +$T$ is in K. Defaults to 0.0. + +`scaling` +: Dimensionless multiplicative scaling factor applied to the rate coefficient. +Defaults to 1.0. + +The parameters `n`, `m`, `E`, and `z` are not permitted for this model. + +This model is intended for detailed V-T and V-V relaxation reactions involving +vibrationally resolved species. All vibrational species in the reaction must belong +to the same vibrational family. + +Example: + +```yaml +- equation: N2(v2) + O => N2(v1) + O + type: vibrational-relaxation + vibration-model: multi-state-resolved + rate-constant: + A: 6.02e+23 cm^3/mol/s + b: 1.0 + B: -34.03 + C: 33.11 + D: 0.0 + scaling: 2.0 +``` + +(sec-vibrational-relaxation-castela-rate)= + +### Vibrational relaxation: castela + +The `castela` vibrational relaxation rate expression is based on the relaxation time + +$$ +\tau_k = +\frac{p_0}{p_k} +\exp \left[ +a_k \left(T^{-1/3} - b_k \right) - 18.42 +\right] +$$ + +and is converted internally to the equivalent bimolecular rate coefficient + +$$ +k_k(T) = +\frac{R T}{p_0} +\exp \left[ +18.42 + a_k b_k - a_k T^{-1/3} +\right] +$$ + +where $T$ is the gas temperature in K, $R$ is the gas constant, and $p_0$ is the +reference pressure. + +The rate parameters are specified as a mapping with fields: + +`a` +: Castela coefficient $a_k$. This value is interpreted assuming that $T$ is in K. + +`b` +: Castela coefficient $b_k$. This value is interpreted assuming that $T$ is in K. + +`reference-pressure` +: Reference pressure $p_0$. Defaults to 1 atm. + +The parameter `A` is not specified by the user for this model. The internal +pre-exponential factor is computed as $R / p_0$. The parameters `A`, `n`, `K`, `B`, +`C`, `D`, `m`, `E`, `z`, and `scaling` are not permitted. + +The `castela` model is only valid for N2 vibrational relaxation and only supports +the colliders `N2`, `O2`, and `O`. It describes relaxation to the ground state, so +the reaction must have the form `N2(v) + M => N2 + M`, where `M` is one of the +supported colliders. + +This model is intended to implement the mean vibrational energy equation model by +the means of a fictitious species N2(v) lumping together all N2 vibrational states. + +Example: + +```yaml +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: castela + rate-constant: + a: 72.4 + b: 0.015 + reference-pressure: 1 atm +``` + +(sec-vibrational-relaxation-starikovskiy-rate)= + +### Vibrational relaxation: starikovskiy + +The `starikovskiy` vibrational relaxation rate expression is + +$$ +k_f = +A T^n +\exp \left( +K ++ B T^{-1/3} ++ C T^{-m} ++ D T^{-z} +\right) +$$ + +where $T$ is the gas temperature in K. + +The rate parameters are specified as a mapping with fields: + +`A` +: The pre-exponential factor. The units are the standard Cantera rate coefficient +units for the reaction order and are converted using the unit system. + +`n` +: The dimensionless temperature exponent. Defaults to 0.0. + +`K` +: Dimensionless constant term in the exponential. Defaults to 0.0. + +`B` +: Signed coefficient multiplying $T^{-1/3}$. This value is interpreted assuming +that $T$ is in K. Defaults to 0.0. + +`C` +: Signed coefficient multiplying $T^{-m}$. This value is interpreted assuming +that $T$ is in K. Defaults to 0.0. + +`m` +: Positive exponent used by the `C` term. Defaults to 1.0. + +`D` +: Signed coefficient multiplying $T^{-z}$. This value is interpreted assuming +that $T$ is in K. Defaults to 0.0. + +`z` +: Positive exponent used by the `D` term. Defaults to 1.0. + +The parameters `b` and `scaling` are not permitted for this model. + +The coefficients `B`, `C`, and `D` are signed values. For example, a negative +coefficient multiplying $T^{-1/3}$ should be written directly as `B: -...` in +the YAML input. + +This model describes relaxation of one vibrationally excited reactant to its ground +state by collision with an unchanged collider: + +```yaml +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: starikovskiy + rate-constant: + A: 6.02214076e+23 cm^3/mol/s + n: 1.0 + K: -34.03 + B: -33.11 + C: 0.0 + m: 1.0 + D: 0.0 + z: 1.0 +``` + +This model is intended to extend the mean vibrational energy equation model from castela +by allowing collisions with more potential colliders. It should be employend alongside a +fictitious species X(v) lumping together all X molecule vibrational states for each molecule +that the user wishes to descibe vibrationally. (sec-yaml-efficiencies)= ## Efficiencies @@ -285,6 +520,51 @@ Example: rate-constant: {A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol} ``` +(sec-yaml-vibrational-relaxation)= +### `vibrational-relaxation` +In plasma phases, vibrationally excited species play an important role in +the process of slow gas heating. Their vibrational-translational relaxation +rate is therefore important to model. This class offers four options to +model VT-relaxation rates: +- [`constant`](sec-vibrational-relaxation-constant-rate) +- [`multi-state-resolved`](sec-vibrational-relaxation-multi-state-resolved-rate) +- [`castela`](sec-vibrational-relaxation-castela-rate) +- [`starikovskiy`](sec-vibrational-relaxation-starikovskiy-rate) + +Example: + +```yaml +- equation: N2(v2) + O => N2(v1) + O + type: vibrational-relaxation + vibration-model: multi-state-resolved + rate-constant: + A: 6.02e+23 + b: 1.0 + B: -34.03 + C: 33.11 + D: 0.0 + scaling: 2.0 + +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: constant + rate-constant: {A: 1e2} + +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: castela + rate-constant: + a: 72.4 + b: 0.015 + reference-pressure: 1 atm + +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: starikovskiy + rate-constant: {A: 6.0221407600e+23, n: 1.0, K: -34.03, B: -33.11, C: 0.0, m: 1.0, D: 0.0, z: 1.0} + +``` + (sec-yaml-electron-collision-plasma)= ### `electron-collision-plasma` From b50eb2fbfd732cab4f36286f369553a5a8868250 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Mon, 13 Jul 2026 08:47:09 +0200 Subject: [PATCH 26/27] correct mispelling in code comments --- src/thermo/PlasmaPhase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index deedd3f8afe..11fd3f77333 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -21,7 +21,7 @@ namespace { const double gamma = sqrt(2 * ElectronCharge / ElectronMass); //! A function to check whether a phase species is actually a vibrational reservoir species - //! If the species is actually a reservoir species, the function sets the basename + //! If the species is actually a reservoir species, the function sets the base name //! at the given baseName adress. bool isVibrationalReservoirName(const string& name, string& baseName) { From 53386455f6667f07c0c6943538e92c4e759b496d Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Mon, 13 Jul 2026 17:37:04 +0200 Subject: [PATCH 27/27] added the T-inv optional parameter in TwoTempPlasmaRate, corrected mispellings and added getParameters function --- doc/sphinx/yaml/reactions.md | 9 +++- include/cantera/kinetics/TwoTempPlasmaRate.h | 56 ++++++++++++++++---- src/kinetics/TwoTempPlasmaRate.cpp | 39 ++++++++++++-- 3 files changed, 89 insertions(+), 15 deletions(-) diff --git a/doc/sphinx/yaml/reactions.md b/doc/sphinx/yaml/reactions.md index 84b51686572..44f80ee7b76 100644 --- a/doc/sphinx/yaml/reactions.md +++ b/doc/sphinx/yaml/reactions.md @@ -129,11 +129,16 @@ a mapping with fields: `Ea-electron` : The activation energy term $E_{a,e}$ that is related to the electron temperature +or a corresponding four-element list. -or a corresponding four-element list. The following are equivalent: +`T-inv` +: Optional gas temperature scale for the term $\exp(-T/T_\mathrm{inv})$. +If omitted or set to 0, this term is not included. + +The following are equivalent: ```yaml -{A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol, b-gas: 0} +{A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol, b-gas: 0, T-inv: 0} [17283, -3.1, -5820 J/mol, 1081 J/mol] ``` (sec-vibrational-relaxation-constant-rate)= diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 581e2899841..2449af1efea 100644 --- a/include/cantera/kinetics/TwoTempPlasmaRate.h +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -47,14 +47,19 @@ struct TwoTempPlasmaData : public ReactionData * activation energy for electron is included. * * @f[ - * k_f = A T^{b_g} T_e^b * exp(-E_{a,g}/RT) * exp(E_{a,e}(T_e - T)/(R T T_e)) + * k_f = A T^{b_g} T_e^b + * \exp\left(-\frac{E_{a,g}}{RT}\right) + * \exp\left(\frac{E_{a,e}(T_e - T)}{R T T_e}\right) + * \exp\left(-\frac{T}{T_\mathrm{inv}}\right) * @f] * * where @f$ T_e @f$ is the electron temperature, @f$ E_{a,g} @f$ is the activation * energy for gas, and @f$ E_{a,e} @f$ is the activation energy for electron, see * Kossyi, et al. @cite kossyi1992. - * The optional gas temperature exponent b_g defaults to zero, which stricly corresponds to @cite kossyi1992. - * If b_g is non-zero, a generalisation is used. + * The optional gas temperature exponent b_g defaults to zero, which strictly corresponds to @cite kossyi1992. + * If b_g is non-zero, a generalisation is used. The optional temperature scale + * @f$ T_\mathrm{inv} @f$ adds the term @f$ \exp(-T/T_\mathrm{inv}) @f$ @cite capitelli2013. + * If @f$ T_\mathrm{inv} = 0 @f$, this term is omitted. * * @ingroup arrheniusGroup */ @@ -65,7 +70,7 @@ class TwoTempPlasmaRate : public ArrheniusBase //! Constructor. /*! - * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units. + * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units * depend on the reaction order and the dimensionality (surface or bulk). * @param b Electron temperature exponent (non-dimensional). * @param Ea Activation energy in energy units [J/kmol]. @@ -77,7 +82,7 @@ class TwoTempPlasmaRate : public ArrheniusBase //! Constructor. /*! - * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units. + * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units * depend on the reaction order and the dimensionality (surface or bulk). * @param b Electron temperature exponent (non-dimensional). * @param Ea Activation energy in energy units [J/kmol]. Defaults to 0. @@ -85,6 +90,21 @@ class TwoTempPlasmaRate : public ArrheniusBase */ TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); + //! Constructor in the cases where T_inv is necessary. + /*! + * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units + * depend on the reaction order and the dimensionality (surface or bulk). + * @param b Electron temperature exponent (non-dimensional). + * @param Ea Activation energy in energy units [J/kmol]. + * @param EE Activation electron energy in energy units [J/kmol]. Defaults to 0. + * @param bg Gas temperature exponent (non-dimensional). Defaults to 0. + * @param Tinv Temperature scale for the term @f$ \exp(-T/T_\mathrm{inv}) @f$ [K]. + * If zero, this term is omitted. Defaults to 0. + * @since New in %Cantera 4.0 + */ + TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg, + double Tinv); + //! Constructor based on an AnyMap object instead of all parameters directly. TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}); @@ -106,10 +126,18 @@ class TwoTempPlasmaRate : public ArrheniusBase * @param shared_data data shared by all reactions of a given type. */ double evalFromStruct(const TwoTempPlasmaData& shared_data) const { - // m_E4_R is the electron activation (in temperature units) - return m_A * std::exp(m_bg * shared_data.logT + m_b * shared_data.logTe - - m_Ea_R * shared_data.recipT + m_E4_R * (shared_data.electronTemp - shared_data.temperature) - * shared_data.recipTe * shared_data.recipT); + // m_E4_R is the electron activation energy in temperature units. + double logRate = m_bg * shared_data.logT + + m_b * shared_data.logTe + - m_Ea_R * shared_data.recipT + + m_E4_R * (shared_data.electronTemp - shared_data.temperature) + * shared_data.recipTe * shared_data.recipT; + + if (m_Tinv != 0.0) { + logRate += -shared_data.temperature / m_Tinv; + } + + return m_A * std::exp(logRate); } //! Evaluate derivative of reaction rate with respect to temperature @@ -127,9 +155,19 @@ class TwoTempPlasmaRate : public ArrheniusBase } protected: + //! Get parameters. + void getParameters(AnyMap& node) const override; + //! Gas temperature exponent. //! @since New in %Cantera 4.0 double m_bg = 0.0; + + //! Temperature scale for the optional term @f$ \exp(-T/T_\mathrm{inv}) @f$. + /*! + * A value of zero disables this term. + * @since New in %Cantera 4.0 + */ + double m_Tinv = 0.0; }; } diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index d97bc4cd88c..a53ad3256d5 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -58,13 +58,20 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) m_E4_str = "Ea-electron"; m_E4_R = EE / GasConstant; m_bg = 0.0; + m_Tinv = 0.0; } TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg) - : ArrheniusBase(A, b, Ea) + : TwoTempPlasmaRate(A, b, Ea, EE, bg, 0.0) +{ +} + +TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE, + double bg, double Tinv) + : TwoTempPlasmaRate(A, b, Ea, EE) { - TwoTempPlasmaRate(A, b, Ea, EE); m_bg = bg; + m_Tinv = Tinv; } TwoTempPlasmaRate::TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units) @@ -72,14 +79,38 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_u { setParameters(node, rate_units); m_bg = node.getDouble("b-gas", 0.0); + m_Tinv = node.getDouble("T-inv", 0.0); +} + +void TwoTempPlasmaRate::getParameters(AnyMap& node) const +{ + ArrheniusBase::getParameters(node); + + auto& rateNode = node["rate-constant"].as(); + + if (m_bg != 0.0) { + rateNode["b-gas"] = m_bg; + } + + if (m_Tinv != 0.0) { + rateNode["T-inv"] = m_Tinv; + } + + rateNode.setFlowStyle(); } double TwoTempPlasmaRate::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const { warn_user("TwoTempPlasmaRate::ddTScaledFromStruct", "Temperature derivative does not consider changes of electron temperature."); - return m_bg * shared_data.recipT - + (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; + double derivative = m_bg * shared_data.recipT + (m_Ea_R - m_E4_R) + * shared_data.recipT * shared_data.recipT; + + if (m_Tinv != 0.0) { + derivative += -1.0 / m_Tinv; + } + + return derivative; } void TwoTempPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin)