From 880c24eec16396f4440290298112f0c445c784be Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 6 Feb 2022 08:41:37 -0600 Subject: [PATCH 01/10] [Kinetics] Introduce Arrhenius3 --- include/cantera/kinetics/Arrhenius.h | 62 +++++++++++++++++++++------- 1 file changed, 46 insertions(+), 16 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 81c4c102d8..4ef9e0e1b0 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -30,11 +30,6 @@ class AnyMap; * This base class provides a minimally functional interface that allows for parameter * access from derived classes as well as classes that use Arrhenius-type expressions * internally, for example FalloffRate and PlogRate. - * - * @todo supersedes Arrhenius2 and will replace Arrhenius(2) after Cantera 2.6, - * The class should be renamed to Arrhenius after removal of Arrhenius2. The new - * behavior can be forced in self-compiled Cantera installations by defining - * CT_NO_LEGACY_REACTIONS_26 via the 'no_legacy_reactions' option in SCons. */ class ArrheniusBase { @@ -165,13 +160,54 @@ class ArrheniusBase * \f] * * @ingroup arrheniusGroup + * + * @todo supersedes Arrhenius2 and will replace Arrhenius(2) after Cantera 2.6, + * The class should be renamed to Arrhenius after removal of Arrhenius2. The new + * behavior can be forced in self-compiled Cantera installations by defining + * CT_NO_LEGACY_REACTIONS_26 via the 'no_legacy_reactions' option in SCons. */ -class ArrheniusRate final : public ArrheniusBase, public ReactionRate +class Arrhenius3 : public ArrheniusBase { public: - ArrheniusRate() = default; using ArrheniusBase::ArrheniusBase; // inherit constructors + const std::string rateType() const { + return "Arrhenius"; + } + + //! Evaluate derivative of reaction rate with respect to temperature + //! divided by reaction rate + /*! + * @internal Non-virtual method that should not be overloaded + */ + double ddTScaled(double logT, double recipT) const { + return (m_Ea_R * recipT + m_b) * recipT; + } + + //! Return the activation energy divided by the gas constant (i.e. the + //! activation temperature) [K] + double activationEnergy_R() const { + return m_Ea_R; + } +}; + + +//! Arrhenius reaction rate type depends only on temperature +/*! + * A reaction rate coefficient of the following form. + * + * \f[ + * k_f = A T^b \exp (-Ea/RT) + * \f] + * + * @ingroup arrheniusGroup + */ +class ArrheniusRate final : public Arrhenius3, public ReactionRate +{ +public: + ArrheniusRate() = default; + using Arrhenius3::Arrhenius3; // inherit constructors + //! Constructor based on AnyMap content ArrheniusRate(const AnyMap& node, const UnitStack& rate_units={}) { setParameters(node, rate_units); @@ -183,7 +219,7 @@ class ArrheniusRate final : public ArrheniusBase, public ReactionRate //! Identifier of reaction rate type virtual const std::string type() const override { - return "Arrhenius"; + return Arrhenius3::rateType(); } //! Perform object setup based on AnyMap node information @@ -204,7 +240,7 @@ class ArrheniusRate final : public ArrheniusBase, public ReactionRate * @param shared_data data shared by all reactions of a given type */ double evalFromStruct(const ReactionData& shared_data) const { - return m_A * std::exp(m_b * shared_data.logT - m_Ea_R * shared_data.recipT); + return evalRate(shared_data.logT, shared_data.recipT); } //! Evaluate derivative of reaction rate with respect to temperature @@ -213,13 +249,7 @@ class ArrheniusRate final : public ArrheniusBase, public ReactionRate * @param shared_data data shared by all reactions of a given type */ virtual double ddTScaledFromStruct(const ReactionData& shared_data) const { - return (m_Ea_R * shared_data.recipT + m_b) * shared_data.recipT; - } - - //! Return the activation energy divided by the gas constant (i.e. the - //! activation temperature) [K] - double activationEnergy_R() const { - return m_Ea_R; + return ddTScaled(shared_data.logT, shared_data.recipT); } }; From bde21edbba5b43b3bec0e2edd79c568854daea21 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 6 Feb 2022 09:28:08 -0600 Subject: [PATCH 02/10] Replace ArrheniusBase by Arrhenius3 as appropriate --- include/cantera/kinetics/Arrhenius.h | 37 +++++++++++--------------- include/cantera/kinetics/Falloff.h | 20 +++++++------- include/cantera/kinetics/RxnRates.h | 23 +++++++--------- interfaces/cython/cantera/_cantera.pxd | 26 +++++++++--------- interfaces/cython/cantera/reaction.pyx | 26 +++++++++--------- src/kinetics/Falloff.cpp | 12 ++++----- src/kinetics/RxnRates.cpp | 32 +++++++++++----------- test/kinetics/kineticsFromScratch3.cpp | 20 +++++++------- 8 files changed, 94 insertions(+), 102 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 4ef9e0e1b0..ad77393ae3 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -69,22 +69,6 @@ class ArrheniusBase //! Check rate expression void checkRate(const std::string& equation, const AnyMap& node); - //! Evaluate reaction rate - /*! - * @internal Non-virtual method that should not be overloaded - */ - double evalRate(double logT, double recipT) const { - return m_A * std::exp(m_b * logT - m_Ea_R * recipT); - } - - //! Evaluate natural logarithm of the rate constant. - /*! - * @internal Non-virtual method that should not be overloaded - */ - double evalLog(double logT, double recipT) const { - return m_logA + m_b * logT - m_Ea_R * recipT; - } - //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending //! on the reaction order) double preExponentialFactor() const { @@ -175,6 +159,22 @@ class Arrhenius3 : public ArrheniusBase return "Arrhenius"; } + //! Evaluate reaction rate + /*! + * @internal Non-virtual method that should not be overloaded + */ + double evalRate(double logT, double recipT) const { + return m_A * std::exp(m_b * logT - m_Ea_R * recipT); + } + + //! Evaluate natural logarithm of the rate constant. + /*! + * @internal Non-virtual method that should not be overloaded + */ + double evalLog(double logT, double recipT) const { + return m_logA + m_b * logT - m_Ea_R * recipT; + } + //! Evaluate derivative of reaction rate with respect to temperature //! divided by reaction rate /*! @@ -222,11 +222,6 @@ class ArrheniusRate final : public Arrhenius3, public ReactionRate return Arrhenius3::rateType(); } - //! Perform object setup based on AnyMap node information - /*! - * @param node AnyMap containing rate information - * @param rate_units Unit definitions specific to rate information - */ virtual void setParameters(const AnyMap& node, const UnitStack& rate_units) override; virtual void getParameters(AnyMap& node) const override; diff --git a/include/cantera/kinetics/Falloff.h b/include/cantera/kinetics/Falloff.h index c18452a7f0..799512cbe8 100644 --- a/include/cantera/kinetics/Falloff.h +++ b/include/cantera/kinetics/Falloff.h @@ -188,24 +188,24 @@ class FalloffRate : public ReactionRate } //! Get reaction rate in the low-pressure limit - ArrheniusBase& lowRate() { + Arrhenius3& lowRate() { return m_lowRate; } //! Set reaction rate in the low-pressure limit - void setLowRate(const ArrheniusBase& low); + void setLowRate(const Arrhenius3& low); //! Get reaction rate in the high-pressure limit - ArrheniusBase& highRate() { + Arrhenius3& highRate() { return m_highRate; } //! Set reaction rate in the high-pressure limit - void setHighRate(const ArrheniusBase& high); + void setHighRate(const Arrhenius3& high); protected: - ArrheniusBase m_lowRate; //!< The reaction rate in the low-pressure limit - ArrheniusBase m_highRate; //!< The reaction rate in the high-pressure limit + Arrhenius3 m_lowRate; //!< The reaction rate in the low-pressure limit + Arrhenius3 m_highRate; //!< The reaction rate in the high-pressure limit bool m_chemicallyActivated; //!< Flag labeling reaction as chemically activated bool m_negativeA_ok; //!< Flag indicating whether negative A values are permitted @@ -234,7 +234,7 @@ class LindemannRate final : public FalloffRate } LindemannRate( - const ArrheniusBase& low, const ArrheniusBase& high, const vector_fp& c) + const Arrhenius3& low, const Arrhenius3& high, const vector_fp& c) : LindemannRate() { m_lowRate = low; @@ -295,7 +295,7 @@ class TroeRate final : public FalloffRate setParameters(node, rate_units); } - TroeRate(const ArrheniusBase& low, const ArrheniusBase& high, const vector_fp& c) + TroeRate(const Arrhenius3& low, const Arrhenius3& high, const vector_fp& c) : TroeRate() { m_lowRate = low; @@ -397,7 +397,7 @@ class SriRate final : public FalloffRate setParameters(node, rate_units); } - SriRate(const ArrheniusBase& low, const ArrheniusBase& high, const vector_fp& c) + SriRate(const Arrhenius3& low, const Arrhenius3& high, const vector_fp& c) : SriRate() { m_lowRate = low; @@ -507,7 +507,7 @@ class TsangRate final : public FalloffRate setParameters(node, rate_units); } - TsangRate(const ArrheniusBase& low, const ArrheniusBase& high, const vector_fp& c) + TsangRate(const Arrhenius3& low, const Arrhenius3& high, const vector_fp& c) : TsangRate() { m_lowRate = low; diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 7e48737513..4d530120b7 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -36,7 +36,7 @@ class Func1; * k_f = A T^b \exp (-E/RT) * \f] */ -class Arrhenius2 : public ArrheniusBase +class Arrhenius2 : public Arrhenius3 { public: //! Default constructor. @@ -57,11 +57,12 @@ class Arrhenius2 : public ArrheniusBase Arrhenius2(const AnyValue& rate, const UnitSystem& units, const Units& rate_units); - Arrhenius2(const ArrheniusBase& other); + //! Converting constructor (to facilitate back-ward compatibility) + Arrhenius2(const Arrhenius3& other); void setRateParameters(const AnyValue& rate, const UnitSystem& units, const Units& rate_units); - using ArrheniusBase::setRateParameters; + using Arrhenius3::setRateParameters; //! Return parameters - two-parameter version void getParameters(AnyMap& node, const Units& rate_units) const; @@ -90,12 +91,6 @@ class Arrhenius2 : public ArrheniusBase doublereal updateRC(doublereal logT, doublereal recipT) const { return m_A * std::exp(m_b*logT - m_Ea_R*recipT); } - - //! Return the activation energy divided by the gas constant (i.e. the - //! activation temperature) [K] - doublereal activationEnergy_R() const { - return m_Ea_R; - } }; @@ -192,7 +187,7 @@ class SurfaceArrhenius #ifdef CT_NO_LEGACY_REACTIONS_26 -typedef ArrheniusBase Arrhenius; +typedef Arrhenius3 Arrhenius; #else typedef Arrhenius2 Arrhenius; #endif @@ -223,7 +218,7 @@ class PlogRate final : public ReactionRate PlogRate(); //! Constructor from Arrhenius rate expressions at a set of pressures - explicit PlogRate(const std::multimap& rates); + explicit PlogRate(const std::multimap& rates); //! Constructor using legacy Arrhenius2 framework explicit PlogRate(const std::multimap& rates); @@ -285,7 +280,7 @@ class PlogRate final : public ReactionRate void setup(const std::multimap& rates); //! Set up Plog object - void setRates(const std::multimap& rates); + void setRates(const std::multimap& rates); //! Update concentration-dependent parts of the rate coefficient. //! @param c natural log of the pressure in Pa @@ -364,14 +359,14 @@ class PlogRate final : public ReactionRate //! Return the pressures and Arrhenius expressions which comprise this //! reaction. - std::multimap getRates() const; + std::multimap getRates() const; protected: //! log(p) to (index range) in the rates_ vector std::map> pressures_; // Rate expressions which are referenced by the indices stored in pressures_ - std::vector rates_; + std::vector rates_; double logP_; //!< log(p) at the current state double logP1_, logP2_; //!< log(p) at the lower / upper pressure reference diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 9340be59eb..1895268fb7 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -459,15 +459,17 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxArrheniusBase "Cantera::ArrheniusBase": CxxArrheniusBase() - CxxArrheniusBase(double, double, double) - double evalRate(double, double) double preExponentialFactor() double temperatureExponent() double activationEnergy() cbool allowNegativePreExponentialFactor() void setAllowNegativePreExponentialFactor(bool) - cdef cppclass CxxArrhenius2 "Cantera::Arrhenius2" (CxxArrheniusBase): + cdef cppclass CxxArrhenius "Cantera::Arrhenius3" (CxxArrheniusBase): + CxxArrhenius(double, double, double) + double evalRate(double, double) + + cdef cppclass CxxArrhenius2 "Cantera::Arrhenius2" (CxxArrhenius): CxxArrhenius2(double, double, double) cdef cppclass CxxReactionRate "Cantera::ReactionRate": @@ -477,7 +479,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": double eval(double, double) except +translate_exception CxxAnyMap parameters() except +translate_exception - cdef cppclass CxxArrheniusRate "Cantera::ArrheniusRate" (CxxReactionRate, CxxArrheniusBase): + cdef cppclass CxxArrheniusRate "Cantera::ArrheniusRate" (CxxReactionRate, CxxArrhenius): CxxArrheniusRate(CxxAnyMap) except +translate_exception CxxArrheniusRate(double, double, double) @@ -499,10 +501,10 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": void setAllowNegativePreExponentialFactor(bool) cbool chemicallyActivated() void setChemicallyActivated(bool) - CxxArrheniusBase& lowRate() - void setLowRate(CxxArrheniusBase&) except +translate_exception - CxxArrheniusBase& highRate() - void setHighRate(CxxArrheniusBase&) except +translate_exception + CxxArrhenius& lowRate() + void setLowRate(CxxArrhenius&) except +translate_exception + CxxArrhenius& highRate() + void setHighRate(CxxArrhenius&) except +translate_exception void getFalloffCoeffs(vector[double]&) void setFalloffCoeffs(vector[double]&) except +translate_exception double evalF(double, double) except +translate_exception @@ -522,8 +524,8 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxPlogRate "Cantera::PlogRate" (CxxReactionRate): CxxPlogRate() CxxPlogRate(CxxAnyMap) except +translate_exception - CxxPlogRate(multimap[double, CxxArrheniusBase]) - multimap[double, CxxArrheniusBase] getRates() + CxxPlogRate(multimap[double, CxxArrhenius]) + multimap[double, CxxArrhenius] getRates() cdef cppclass CxxChebyshevRate "Cantera::ChebyshevRate" (CxxReactionRate): CxxChebyshevRate() @@ -1404,11 +1406,11 @@ cdef class CustomReaction(Reaction): cdef class Arrhenius: cdef CxxArrhenius2* legacy # used by legacy objects only - cdef CxxArrheniusBase* base + cdef CxxArrhenius* base cdef cbool own_rate cdef Reaction reaction # parent reaction, to prevent garbage collection @staticmethod - cdef wrap(CxxArrheniusBase*) + cdef wrap(CxxArrhenius*) cdef class BlowersMasel: cdef CxxBlowersMasel2* rate diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index b967dc2417..22190d73e6 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -214,8 +214,8 @@ cdef class ArrheniusRate(ArrheniusTypeRate): cdef set_cxx_object(self): self.rate = self._rate.get() - # CxxArrheniusBase does not have a common base with CxxReactionRate - self.base = self.cxx_object() + # CxxArrhenius does not have a common base with CxxReactionRate + self.base = self.cxx_object() cdef CxxArrheniusRate* cxx_object(self): return self.rate @@ -522,17 +522,17 @@ cdef class PlogRate(ReactionRate): """ def __get__(self): rates = [] - cdef multimap[double, CxxArrheniusBase] cxxrates - cdef pair[double, CxxArrheniusBase] p_rate + cdef multimap[double, CxxArrhenius] cxxrates + cdef pair[double, CxxArrhenius] p_rate cxxrates = self.cxx_object().getRates() for p_rate in cxxrates: - rates.append((p_rate.first, copyArrheniusBase(&p_rate.second))) + rates.append((p_rate.first, copyArrhenius(&p_rate.second))) return rates def __set__(self, rates): - cdef multimap[double, CxxArrheniusBase] ratemap + cdef multimap[double, CxxArrhenius] ratemap cdef Arrhenius rate - cdef pair[double, CxxArrheniusBase] item + cdef pair[double, CxxArrhenius] item for p, rate in rates: item.first = p item.second = deref(rate.base) @@ -1453,7 +1453,7 @@ cdef class Arrhenius: """ def __cinit__(self, A=0, b=0, E=0, init=True): if init: - self.base = new CxxArrheniusBase(A, b, E) + self.base = new CxxArrhenius(A, b, E) self.own_rate = True self.reaction = None else: @@ -1464,7 +1464,7 @@ cdef class Arrhenius: del self.base @staticmethod - cdef wrap(CxxArrheniusBase* rate): + cdef wrap(CxxArrhenius* rate): r = Arrhenius(init=False) r.base = rate r.reaction = None @@ -1501,7 +1501,7 @@ cdef class Arrhenius: return self.base.evalRate(np.log(T), 1/T) -cdef wrapArrhenius(CxxArrheniusBase* rate, Reaction reaction): +cdef wrapArrhenius(CxxArrhenius* rate, Reaction reaction): r = Arrhenius(init=False) r.base = rate if reaction.uses_legacy: @@ -1510,12 +1510,12 @@ cdef wrapArrhenius(CxxArrheniusBase* rate, Reaction reaction): r.reaction = reaction return r -cdef copyArrhenius(CxxArrhenius2* rate): +cdef copyArrhenius2(CxxArrhenius2* rate): r = Arrhenius(rate.preExponentialFactor(), rate.temperatureExponent(), rate.activationEnergy()) return r -cdef copyArrheniusBase(CxxArrheniusBase* rate): +cdef copyArrhenius(CxxArrhenius* rate): r = Arrhenius(rate.preExponentialFactor(), rate.temperatureExponent(), rate.activationEnergy()) return r @@ -2118,7 +2118,7 @@ cdef class PlogReaction(Reaction): cdef pair[double,CxxArrhenius2] p_rate rates = [] for p_rate in cxxrates: - rates.append((p_rate.first,copyArrhenius(&p_rate.second))) + rates.append((p_rate.first,copyArrhenius2(&p_rate.second))) return rates cdef _legacy_set_rates(self, list rates): diff --git a/src/kinetics/Falloff.cpp b/src/kinetics/Falloff.cpp index e0cdbc77d9..ec5d18ba53 100644 --- a/src/kinetics/Falloff.cpp +++ b/src/kinetics/Falloff.cpp @@ -21,9 +21,9 @@ void FalloffRate::init(const vector_fp& c) setFalloffCoeffs(c); } -void FalloffRate::setLowRate(const ArrheniusBase& low) +void FalloffRate::setLowRate(const Arrhenius3& low) { - ArrheniusBase _low = low; + Arrhenius3 _low = low; _low.setAllowNegativePreExponentialFactor(m_negativeA_ok); _low.checkRate("", AnyMap()); if (_low.preExponentialFactor() * m_highRate.preExponentialFactor() < 0.) { @@ -34,9 +34,9 @@ void FalloffRate::setLowRate(const ArrheniusBase& low) m_lowRate = std::move(_low); } -void FalloffRate::setHighRate(const ArrheniusBase& high) +void FalloffRate::setHighRate(const Arrhenius3& high) { - ArrheniusBase _high = high; + Arrhenius3 _high = high; _high.setAllowNegativePreExponentialFactor(m_negativeA_ok); _high.checkRate("", AnyMap()); if (m_lowRate.preExponentialFactor() * _high.preExponentialFactor() < 0.) { @@ -83,12 +83,12 @@ void FalloffRate::setParameters(const AnyMap& node, const UnitStack& rate_units) } } if (node.hasKey("low-P-rate-constant")) { - m_lowRate = ArrheniusBase( + m_lowRate = Arrhenius3( node["low-P-rate-constant"], node.units(), low_rate_units); m_lowRate.setAllowNegativePreExponentialFactor(m_negativeA_ok); } if (node.hasKey("high-P-rate-constant")) { - m_highRate = ArrheniusBase( + m_highRate = Arrhenius3( node["high-P-rate-constant"], node.units(), high_rate_units); m_highRate.setAllowNegativePreExponentialFactor(m_negativeA_ok); } diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 6908431f03..948da372fc 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -9,7 +9,7 @@ namespace Cantera { Arrhenius2::Arrhenius2() - : ArrheniusBase() + : Arrhenius3() { m_b = 0.0; m_A = 0.0; @@ -17,7 +17,7 @@ Arrhenius2::Arrhenius2() } Arrhenius2::Arrhenius2(doublereal A, doublereal b, doublereal E) - : ArrheniusBase(A, b, E * GasConstant) + : Arrhenius3(A, b, E * GasConstant) { if (m_A <= 0.0) { m_logA = -1.0E300; @@ -30,10 +30,10 @@ Arrhenius2::Arrhenius2(const AnyValue& rate, setRateParameters(rate, units, rate_units); } -Arrhenius2::Arrhenius2(const ArrheniusBase& other) - : ArrheniusBase(other.preExponentialFactor(), - other.temperatureExponent(), - other.activationEnergy()) +Arrhenius2::Arrhenius2(const Arrhenius3& other) + : Arrhenius3(other.preExponentialFactor(), + other.temperatureExponent(), + other.activationEnergy()) { } @@ -41,7 +41,7 @@ void Arrhenius2::setRateParameters(const AnyValue& rate, const UnitSystem& units, const Units& rate_units) { UnitStack units_stack(rate_units); - ArrheniusBase::setRateParameters(rate, units, units_stack); + Arrhenius3::setRateParameters(rate, units, units_stack); if (m_A <= 0.0) { m_logA = -1.0E300; } @@ -103,7 +103,7 @@ PlogRate::PlogRate() { } -PlogRate::PlogRate(const std::multimap& rates) +PlogRate::PlogRate(const std::multimap& rates) : PlogRate() { setRates(rates); @@ -130,16 +130,16 @@ void PlogRate::setParameters(const AnyMap& node, const UnitStack& units) void PlogRate::setRateParameters(const std::vector& rates, const UnitSystem& units, const Units& rate_units) { - std::multimap multi_rates; + std::multimap multi_rates; if (rates.size()) { for (const auto& rate : rates) { multi_rates.insert({rate.convert("P", "Pa"), - ArrheniusBase(AnyValue(rate), units, rate_units)}); + Arrhenius3(AnyValue(rate), units, rate_units)}); } } else { // ensure that reaction rate can be evaluated (but returns NaN) - multi_rates.insert({1.e-7, ArrheniusBase(NAN, NAN, NAN)}); - multi_rates.insert({1.e14, ArrheniusBase(NAN, NAN, NAN)}); + multi_rates.insert({1.e-7, Arrhenius3(NAN, NAN, NAN)}); + multi_rates.insert({1.e14, Arrhenius3(NAN, NAN, NAN)}); } setRates(multi_rates); } @@ -170,14 +170,14 @@ void PlogRate::getParameters(AnyMap& rateNode, const Units& rate_units) const void PlogRate::setup(const std::multimap& rates) { - std::multimap rates2; + std::multimap rates2; for (const auto& item : rates) { rates2.emplace(item.first, item.second); } setRates(rates2); } -void PlogRate::setRates(const std::multimap& rates) +void PlogRate::setRates(const std::multimap& rates) { size_t j = 0; rates_.clear(); @@ -238,9 +238,9 @@ std::vector > PlogRate::rates() const return out; } -std::multimap PlogRate::getRates() const +std::multimap PlogRate::getRates() const { - std::multimap rateMap; + std::multimap rateMap; // initial preincrement to skip rate for P --> 0 for (auto iter = ++pressures_.begin(); iter->first < 1000; // skip rates for (P --> infinity) diff --git a/test/kinetics/kineticsFromScratch3.cpp b/test/kinetics/kineticsFromScratch3.cpp index 51fa0c5c80..6c4132c64c 100644 --- a/test/kinetics/kineticsFromScratch3.cpp +++ b/test/kinetics/kineticsFromScratch3.cpp @@ -136,11 +136,11 @@ TEST_F(KineticsFromScratch3, add_plog_reaction) // [(100.0, 'atm'), 5.963200e+56, -11.529, 52599.6]) Composition reac = parseCompString("H2:1, O2:1"); Composition prod = parseCompString("OH:2"); - std::multimap rates { - { 0.01*101325, ArrheniusBase(1.212400e+16, -0.5779, 10872.7 * 4184.0) }, - { 1.0*101325, ArrheniusBase(4.910800e+31, -4.8507, 24772.8 * 4184.0) }, - { 10.0*101325, ArrheniusBase(1.286600e+47, -9.0246, 39796.5 * 4184.0) }, - { 100.0*101325, ArrheniusBase(5.963200e+56, -11.529, 52599.6 * 4184.0) } + std::multimap rates { + { 0.01*101325, Arrhenius3(1.212400e+16, -0.5779, 10872.7 * 4184.0) }, + { 1.0*101325, Arrhenius3(4.910800e+31, -4.8507, 24772.8 * 4184.0) }, + { 10.0*101325, Arrhenius3(1.286600e+47, -9.0246, 39796.5 * 4184.0) }, + { 100.0*101325, Arrhenius3(5.963200e+56, -11.529, 52599.6 * 4184.0) } }; auto R = make_shared(reac, prod, make_shared(rates)); @@ -152,11 +152,11 @@ TEST_F(KineticsFromScratch3, plog_invalid_rate) { Composition reac = parseCompString("H2:1, O2:1"); Composition prod = parseCompString("OH:2"); - std::multimap rates { - { 0.01*101325, ArrheniusBase(1.2124e+16, -0.5779, 10872.7 * 4184.0) }, - { 10.0*101325, ArrheniusBase(1e15, -1, 10000 * 4184.0) }, - { 10.0*101325, ArrheniusBase(-2e20, -2.0, 20000 * 4184.0) }, - { 100.0*101325, ArrheniusBase(5.9632e+56, -11.529, 52599.6 * 4184.0) } + std::multimap rates { + { 0.01*101325, Arrhenius3(1.2124e+16, -0.5779, 10872.7 * 4184.0) }, + { 10.0*101325, Arrhenius3(1e15, -1, 10000 * 4184.0) }, + { 10.0*101325, Arrhenius3(-2e20, -2.0, 20000 * 4184.0) }, + { 100.0*101325, Arrhenius3(5.9632e+56, -11.529, 52599.6 * 4184.0) } }; auto R = make_shared(reac, prod, make_shared(rates)); From 342981e127038ba6f0a58c0ec4c13e89d6341a38 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 6 Feb 2022 08:03:10 -0600 Subject: [PATCH 03/10] [Kinetics] Re-introduce BlowersMasel class --- include/cantera/kinetics/Arrhenius.h | 125 ++++++++++++++++++--------- src/kinetics/Arrhenius.cpp | 10 +-- 2 files changed, 88 insertions(+), 47 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index ad77393ae3..637f57a3f7 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -372,11 +372,11 @@ class TwoTempPlasmaRate final : public ArrheniusBase, public ReactionRate * * @ingroup arrheniusGroup */ -class BlowersMaselRate final : public ArrheniusBase, public ReactionRate +class BlowersMasel : public ArrheniusBase { public: //! Default constructor. - BlowersMaselRate(); + BlowersMasel(); //! Constructor. /*! @@ -387,46 +387,19 @@ class BlowersMaselRate final : public ArrheniusBase, public ReactionRate * @param w Average bond dissociation energy of the bond being formed and * broken in the reaction, in energy units [J/kmol] */ - BlowersMaselRate(double A, double b, double Ea0, double w); + BlowersMasel(double A, double b, double Ea0, double w); - unique_ptr newMultiRate() const override { - return unique_ptr( - new MultiRate); - } - - //! Constructor based on AnyMap content - BlowersMaselRate(const AnyMap& node, const UnitStack& rate_units={}) - : BlowersMaselRate() - { - setParameters(node, rate_units); - } - - //! Identifier of reaction rate type - virtual const std::string type() const override { + const std::string rateType() const { return "Blowers-Masel"; } - //! Perform object setup based on AnyMap node information - /*! - * @param node AnyMap containing rate information - * @param rate_units Unit definitions specific to rate information - */ - virtual void setParameters( - const AnyMap& node, const UnitStack& rate_units) override; - - virtual void getParameters(AnyMap& node) const override; - - void check(const std::string& equation, const AnyMap& node) override { - checkRate(equation, node); - } - - virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; + virtual void setRateContext(const Reaction& rxn, const Kinetics& kin); //! Update information specific to reaction /*! * @param shared_data data shared by all reactions of a given type */ - void updateFromStruct(const BlowersMaselData& shared_data) { + void updateRate(const BlowersMaselData& shared_data) { if (shared_data.ready) { m_deltaH_R = 0.; for (const auto& item : m_stoich_coeffs) { @@ -440,22 +413,22 @@ class BlowersMaselRate final : public ArrheniusBase, public ReactionRate //! Evaluate reaction rate /*! - * @param shared_data data shared by all reactions of a given type + * @internal Non-virtual method that should not be overloaded */ - double evalFromStruct(const BlowersMaselData& shared_data) { + double evalRate(double logT, double recipT) const { double Ea_R = activationEnergy_R(m_deltaH_R); - return m_A * std::exp(m_b * shared_data.logT - Ea_R * shared_data.recipT); + return m_A * std::exp(m_b * logT - Ea_R * recipT); } //! Evaluate derivative of reaction rate with respect to temperature //! divided by reaction rate /*! - * This method is used to override the numerical derivative, which does not - * consider potential changes due to a changed reaction enthalpy. A corresponding - * warning is raised. - * @param shared_data data shared by all reactions of a given type + * This method does not consider potential changes due to a changed reaction + * enthalpy. A corresponding warning is raised. + * + * @internal Non-virtual method that should not be overloaded */ - virtual double ddTScaledFromStruct(const BlowersMaselData& shared_data) const; + double ddTScaled(double logT, double recipT) const; //! Return the effective activation energy (a function of the delta H of reaction) //! divided by the gas constant (i.e. the activation temperature) [K] @@ -493,7 +466,75 @@ class BlowersMaselRate final : public ArrheniusBase, public ReactionRate double m_deltaH_R; //!< enthalpy change of reaction (in temperature units) }; -} +class BlowersMaselRate final : public BlowersMasel, public ReactionRate +{ +public: + using BlowersMasel::BlowersMasel; // inherit constructors + + unique_ptr newMultiRate() const { + return unique_ptr( + new MultiRate); + } + + //! Constructor based on AnyMap content + BlowersMaselRate(const AnyMap& node, const UnitStack& rate_units={}) + : BlowersMasel() + { + setParameters(node, rate_units); + } + + //! Identifier of reaction rate type + virtual const std::string type() const { + return BlowersMasel::rateType(); + } + + //! Perform object setup based on AnyMap node information + /*! + * @param node AnyMap containing rate information + * @param rate_units Unit definitions specific to rate information + */ + virtual void setParameters(const AnyMap& node, const UnitStack& rate_units); + + virtual void getParameters(AnyMap& node) const; + + void check(const std::string& equation, const AnyMap& node) override { + checkRate(equation, node); + } + + virtual void setContext(const Reaction& rxn, const Kinetics& kin) override { + setRateContext(rxn, kin); + } + + //! Update information specific to reaction + /*! + * @param shared_data data shared by all reactions of a given type + */ + void updateFromStruct(const BlowersMaselData& shared_data) { + updateRate(shared_data); + } + + //! Evaluate reaction rate + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const BlowersMaselData& shared_data) { + return evalRate(shared_data.logT, shared_data.recipT); + } + + //! Evaluate derivative of reaction rate with respect to temperature + //! divided by reaction rate + /*! + * This method is used to override the numerical derivative, which does not + * consider potential changes due to a changed reaction enthalpy. A corresponding + * warning is raised. + * @param shared_data data shared by all reactions of a given type + */ + virtual double ddTScaledFromStruct(const BlowersMaselData& shared_data) const { + return ddTScaled(shared_data.logT, shared_data.recipT); + } +}; + +} #endif diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index c77cdfdec8..cdedb95982 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -189,14 +189,14 @@ void TwoTempPlasmaRate::getParameters(AnyMap& rateNode) const rateNode["type"] = type(); } -BlowersMaselRate::BlowersMaselRate() +BlowersMasel::BlowersMasel() : m_deltaH_R(0.) { m_Ea_str = "Ea0"; m_E4_str = "w"; } -BlowersMaselRate::BlowersMaselRate(double A, double b, double Ea0, double w) +BlowersMasel::BlowersMasel(double A, double b, double Ea0, double w) : ArrheniusBase(A, b, Ea0) , m_deltaH_R(0.) { @@ -205,12 +205,12 @@ BlowersMaselRate::BlowersMaselRate(double A, double b, double Ea0, double w) m_E4_R = w / GasConstant; } -double BlowersMaselRate::ddTScaledFromStruct(const BlowersMaselData& shared_data) const +double BlowersMasel::ddTScaled(double logT, double recipT) const { warn_user("BlowersMaselRate::ddTScaledFromStruct", "Temperature derivative does not consider changes of reaction enthalpy."); double Ea_R = activationEnergy_R(m_deltaH_R); - return m_A * std::exp(m_b * shared_data.logT - Ea_R * shared_data.recipT); + return m_A * std::exp(m_b * logT - Ea_R * recipT); } void BlowersMaselRate::setParameters(const AnyMap& node, const UnitStack& rate_units) @@ -237,7 +237,7 @@ void BlowersMaselRate::getParameters(AnyMap& rateNode) const rateNode["type"] = type(); } -void BlowersMaselRate::setContext(const Reaction& rxn, const Kinetics& kin) +void BlowersMasel::setRateContext(const Reaction& rxn, const Kinetics& kin) { m_stoich_coeffs.clear(); for (const auto& sp : rxn.reactants) { From 3d521aa7325c1bb3bba6049073a101c8353b1bbb Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 6 Feb 2022 10:48:18 -0600 Subject: [PATCH 04/10] [Kinetics] Introduce BulkRate template ArrheniusRate is defined as BulkRate --- include/cantera/kinetics/Arrhenius.h | 147 ++++++++++++++------------- src/kinetics/Arrhenius.cpp | 24 ----- 2 files changed, 79 insertions(+), 92 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 637f57a3f7..1224e30e6d 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -160,28 +160,26 @@ class Arrhenius3 : public ArrheniusBase } //! Evaluate reaction rate - /*! - * @internal Non-virtual method that should not be overloaded - */ + //! @internal Non-virtual method that should not be overloaded double evalRate(double logT, double recipT) const { return m_A * std::exp(m_b * logT - m_Ea_R * recipT); } //! Evaluate natural logarithm of the rate constant. - /*! - * @internal Non-virtual method that should not be overloaded - */ + //! @internal Non-virtual method that should not be overloaded double evalLog(double logT, double recipT) const { return m_logA + m_b * logT - m_Ea_R * recipT; } + //! Evaluate reaction rate + double evalRate(const ArrheniusData& shared_data) const { + return m_A * std::exp(m_b * shared_data.logT - m_Ea_R * shared_data.recipT); + } + //! Evaluate derivative of reaction rate with respect to temperature //! divided by reaction rate - /*! - * @internal Non-virtual method that should not be overloaded - */ - double ddTScaled(double logT, double recipT) const { - return (m_Ea_R * recipT + m_b) * recipT; + double ddTScaled(const ArrheniusData& shared_data) const { + return (m_Ea_R * shared_data.recipT + m_b) * shared_data.recipT; } //! Return the activation energy divided by the gas constant (i.e. the @@ -192,63 +190,6 @@ class Arrhenius3 : public ArrheniusBase }; -//! Arrhenius reaction rate type depends only on temperature -/*! - * A reaction rate coefficient of the following form. - * - * \f[ - * k_f = A T^b \exp (-Ea/RT) - * \f] - * - * @ingroup arrheniusGroup - */ -class ArrheniusRate final : public Arrhenius3, public ReactionRate -{ -public: - ArrheniusRate() = default; - using Arrhenius3::Arrhenius3; // inherit constructors - - //! Constructor based on AnyMap content - ArrheniusRate(const AnyMap& node, const UnitStack& rate_units={}) { - setParameters(node, rate_units); - } - - unique_ptr newMultiRate() const override { - return unique_ptr(new MultiRate); - } - - //! Identifier of reaction rate type - virtual const std::string type() const override { - return Arrhenius3::rateType(); - } - - virtual void setParameters(const AnyMap& node, const UnitStack& rate_units) override; - - virtual void getParameters(AnyMap& node) const override; - - void check(const std::string& equation, const AnyMap& node) override { - checkRate(equation, node); - } - - //! Evaluate reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double evalFromStruct(const ReactionData& shared_data) const { - return evalRate(shared_data.logT, shared_data.recipT); - } - - //! Evaluate derivative of reaction rate with respect to temperature - //! divided by reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - virtual double ddTScaledFromStruct(const ReactionData& shared_data) const { - return ddTScaled(shared_data.logT, shared_data.recipT); - } -}; - - //! Two temperature plasma reaction rate type depends on both //! gas temperature and electron temperature. /*! @@ -535,6 +476,76 @@ class BlowersMaselRate final : public BlowersMasel, public ReactionRate } }; + +//! Template for bulk phase reaction rate specifications +template +class BulkRate final : public RateType, public ReactionRate +{ +public: + BulkRate() = default; + using RateType::RateType; // inherit constructors + + //! Constructor based on AnyMap content + BulkRate(const AnyMap& node, const UnitStack& rate_units={}) { + setParameters(node, rate_units); + } + + unique_ptr newMultiRate() const override { + return unique_ptr( + new MultiRate, DataType>); + } + + virtual const std::string type() const override { + return RateType::rateType(); + } + + virtual void setParameters( + const AnyMap& node, const UnitStack& rate_units) override + { + RateType::m_negativeA_ok = node.getBool("negative-A", false); + if (!node.hasKey("rate-constant")) { + RateType::setRateParameters(AnyValue(), node.units(), rate_units); + return; + } + RateType::setRateParameters(node["rate-constant"], node.units(), rate_units); + } + + virtual void getParameters(AnyMap& node) const override { + if (RateType::m_negativeA_ok) { + node["negative-A"] = true; + } + AnyMap rateNode; + RateType::getRateParameters(rateNode); + if (!rateNode.empty()) { + // RateType object is configured + node["rate-constant"] = std::move(rateNode); + } + } + + void check(const std::string& equation, const AnyMap& node) override { + checkRate(equation, node); + } + + //! Evaluate reaction rate + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const DataType& shared_data) const { + return evalRate(shared_data); + } + + //! Evaluate derivative of reaction rate with respect to temperature + //! divided by reaction rate + /*! + * @param shared_data data shared by all reactions of a given type + */ + double ddTScaledFromStruct(const DataType& shared_data) const { + return ddTScaled(shared_data); + } +}; + +typedef BulkRate ArrheniusRate; + } #endif diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index cdedb95982..6c9b557562 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -126,30 +126,6 @@ void ArrheniusBase::checkRate(const std::string& equation, const AnyMap& node) } } -void ArrheniusRate::setParameters(const AnyMap& node, const UnitStack& rate_units) -{ - m_negativeA_ok = node.getBool("negative-A", false); - if (!node.hasKey("rate-constant")) { - ArrheniusBase::setRateParameters(AnyValue(), node.units(), rate_units); - return; - } - - ArrheniusBase::setRateParameters(node["rate-constant"], node.units(), rate_units); -} - -void ArrheniusRate::getParameters(AnyMap& rateNode) const -{ - if (m_negativeA_ok) { - rateNode["negative-A"] = true; - } - AnyMap node; - ArrheniusBase::getRateParameters(node); - if (!node.empty()) { - // Arrhenius object is configured - rateNode["rate-constant"] = std::move(node); - } -} - TwoTempPlasmaRate::TwoTempPlasmaRate() : ArrheniusBase() { From 7a68a114b5bb9362be11f82842d348feb7e0f3e7 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 6 Feb 2022 11:05:40 -0600 Subject: [PATCH 05/10] [Kinetics] Define BlowersMaselRate as BulkRate --- include/cantera/kinetics/Arrhenius.h | 121 ++++++------------------ include/cantera/kinetics/ReactionRate.h | 4 +- src/kinetics/Arrhenius.cpp | 32 +------ test/kinetics/kineticsFromYaml.cpp | 6 +- 4 files changed, 39 insertions(+), 124 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 1224e30e6d..9eca8673d6 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -88,6 +88,12 @@ class ArrheniusBase return m_Ea_R * GasConstant; } + //! Return the activation energy divided by the gas constant (i.e. the + //! activation temperature) [K] + double activationEnergy_R() const { + return m_Ea_R; + } + // Return units of the reaction rate expression const Units& rateUnits() const { return m_rate_units; @@ -159,6 +165,10 @@ class Arrhenius3 : public ArrheniusBase return "Arrhenius"; } + //! Context (unused) + void setRateContext(const Reaction& rxn, const Kinetics& kin) { + } + //! Evaluate reaction rate //! @internal Non-virtual method that should not be overloaded double evalRate(double logT, double recipT) const { @@ -181,12 +191,6 @@ class Arrhenius3 : public ArrheniusBase double ddTScaled(const ArrheniusData& shared_data) const { return (m_Ea_R * shared_data.recipT + m_b) * shared_data.recipT; } - - //! Return the activation energy divided by the gas constant (i.e. the - //! activation temperature) [K] - double activationEnergy_R() const { - return m_Ea_R; - } }; @@ -334,13 +338,11 @@ class BlowersMasel : public ArrheniusBase return "Blowers-Masel"; } - virtual void setRateContext(const Reaction& rxn, const Kinetics& kin); + void setRateContext(const Reaction& rxn, const Kinetics& kin); //! Update information specific to reaction - /*! - * @param shared_data data shared by all reactions of a given type - */ - void updateRate(const BlowersMaselData& shared_data) { + //! @todo Move to BulkRate (as template function via std::enable_if) + void updateFromStruct(const BlowersMaselData& shared_data) { if (shared_data.ready) { m_deltaH_R = 0.; for (const auto& item : m_stoich_coeffs) { @@ -353,12 +355,10 @@ class BlowersMasel : public ArrheniusBase } //! Evaluate reaction rate - /*! - * @internal Non-virtual method that should not be overloaded - */ - double evalRate(double logT, double recipT) const { - double Ea_R = activationEnergy_R(m_deltaH_R); - return m_A * std::exp(m_b * logT - Ea_R * recipT); + //! @internal Non-virtual method that should not be overloaded + double evalRate(const BlowersMaselData& shared_data) const { + double Ea_R = effectiveActivationEnergy_R(m_deltaH_R); + return m_A * std::exp(m_b * shared_data.logT - Ea_R * shared_data.recipT); } //! Evaluate derivative of reaction rate with respect to temperature @@ -366,14 +366,13 @@ class BlowersMasel : public ArrheniusBase /*! * This method does not consider potential changes due to a changed reaction * enthalpy. A corresponding warning is raised. - * * @internal Non-virtual method that should not be overloaded */ - double ddTScaled(double logT, double recipT) const; + double ddTScaled(const BlowersMaselData& shared_data) const; //! Return the effective activation energy (a function of the delta H of reaction) //! divided by the gas constant (i.e. the activation temperature) [K] - double activationEnergy_R(double deltaH_R) const { + double effectiveActivationEnergy_R(double deltaH_R) const { if (deltaH_R < -4 * m_Ea_R) { return 0.; } @@ -392,7 +391,7 @@ class BlowersMasel : public ArrheniusBase * @param deltaH Enthalpy change of reaction [J/kmol] */ double effectiveActivationEnergy(double deltaH) const { - return activationEnergy_R(deltaH / GasConstant) * GasConstant; + return effectiveActivationEnergy_R(deltaH / GasConstant) * GasConstant; } //! Return the bond dissociation energy *w* [J/kmol] @@ -408,76 +407,7 @@ class BlowersMasel : public ArrheniusBase }; -class BlowersMaselRate final : public BlowersMasel, public ReactionRate -{ -public: - using BlowersMasel::BlowersMasel; // inherit constructors - - unique_ptr newMultiRate() const { - return unique_ptr( - new MultiRate); - } - - //! Constructor based on AnyMap content - BlowersMaselRate(const AnyMap& node, const UnitStack& rate_units={}) - : BlowersMasel() - { - setParameters(node, rate_units); - } - - //! Identifier of reaction rate type - virtual const std::string type() const { - return BlowersMasel::rateType(); - } - - //! Perform object setup based on AnyMap node information - /*! - * @param node AnyMap containing rate information - * @param rate_units Unit definitions specific to rate information - */ - virtual void setParameters(const AnyMap& node, const UnitStack& rate_units); - - virtual void getParameters(AnyMap& node) const; - - void check(const std::string& equation, const AnyMap& node) override { - checkRate(equation, node); - } - - virtual void setContext(const Reaction& rxn, const Kinetics& kin) override { - setRateContext(rxn, kin); - } - - //! Update information specific to reaction - /*! - * @param shared_data data shared by all reactions of a given type - */ - void updateFromStruct(const BlowersMaselData& shared_data) { - updateRate(shared_data); - } - - //! Evaluate reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double evalFromStruct(const BlowersMaselData& shared_data) { - return evalRate(shared_data.logT, shared_data.recipT); - } - - //! Evaluate derivative of reaction rate with respect to temperature - //! divided by reaction rate - /*! - * This method is used to override the numerical derivative, which does not - * consider potential changes due to a changed reaction enthalpy. A corresponding - * warning is raised. - * @param shared_data data shared by all reactions of a given type - */ - virtual double ddTScaledFromStruct(const BlowersMaselData& shared_data) const { - return ddTScaled(shared_data.logT, shared_data.recipT); - } -}; - - -//! Template for bulk phase reaction rate specifications +//! A class template for bulk phase reaction rate specifications template class BulkRate final : public RateType, public ReactionRate { @@ -520,12 +450,20 @@ class BulkRate final : public RateType, public ReactionRate // RateType object is configured node["rate-constant"] = std::move(rateNode); } + if (RateType::rateType() != "Arrhenius") { + node["type"] = type(); + } } void check(const std::string& equation, const AnyMap& node) override { checkRate(equation, node); } + virtual void setContext(const Reaction& rxn, const Kinetics& kin) override { + // as this method is virtual, it cannot be templated + setRateContext(rxn, kin); + } + //! Evaluate reaction rate /*! * @param shared_data data shared by all reactions of a given type @@ -545,6 +483,7 @@ class BulkRate final : public RateType, public ReactionRate }; typedef BulkRate ArrheniusRate; +typedef BulkRate BlowersMaselRate; } diff --git a/include/cantera/kinetics/ReactionRate.h b/include/cantera/kinetics/ReactionRate.h index 19cbba8a57..ca722e5780 100644 --- a/include/cantera/kinetics/ReactionRate.h +++ b/include/cantera/kinetics/ReactionRate.h @@ -117,8 +117,8 @@ class ReactionRate //! Set context of reaction rate evaluation //! @param rxn Reaction object associated with rate //! @param kin Kinetics object used for rate evaluation - //! This method allows for passing of information when a ReactionRate is added - //! to Kinetics a MultiRate reaction evaluator. + //! This method allows for passing of information specific to the associated + //! reaction when a ReactionRate object is added a MultiRate reaction evaluator. virtual void setContext(const Reaction& rxn, const Kinetics& kin) { } diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index 6c9b557562..c6319b6174 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -181,36 +181,12 @@ BlowersMasel::BlowersMasel(double A, double b, double Ea0, double w) m_E4_R = w / GasConstant; } -double BlowersMasel::ddTScaled(double logT, double recipT) const +double BlowersMasel::ddTScaled(const BlowersMaselData& shared_data) const { - warn_user("BlowersMaselRate::ddTScaledFromStruct", + warn_user("BlowersMasel::ddTScaledFromStruct", "Temperature derivative does not consider changes of reaction enthalpy."); - double Ea_R = activationEnergy_R(m_deltaH_R); - return m_A * std::exp(m_b * logT - Ea_R * recipT); -} - -void BlowersMaselRate::setParameters(const AnyMap& node, const UnitStack& rate_units) -{ - m_negativeA_ok = node.getBool("negative-A", false); - if (!node.hasKey("rate-constant")) { - setRateParameters(AnyValue(), node.units(), rate_units); - return; - } - setRateParameters(node["rate-constant"], node.units(), rate_units); -} - -void BlowersMaselRate::getParameters(AnyMap& rateNode) const -{ - if (m_negativeA_ok) { - rateNode["negative-A"] = true; - } - AnyMap node; - ArrheniusBase::getRateParameters(node); - if (!node.empty()) { - // object is configured - rateNode["rate-constant"] = std::move(node); - } - rateNode["type"] = type(); + double Ea_R = effectiveActivationEnergy_R(m_deltaH_R); + return (Ea_R * shared_data.recipT + m_b) * shared_data.recipT; } void BlowersMasel::setRateContext(const Reaction& rxn, const Kinetics& kin) diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 5993db9e07..253a601738 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -289,9 +289,9 @@ TEST(Reaction, BlowersMaselFromYaml) EXPECT_DOUBLE_EQ(rate->preExponentialFactor(), -38.7); EXPECT_DOUBLE_EQ(rate->activationEnergy(), E_intrinsic); EXPECT_DOUBLE_EQ(rate->bondEnergy(), w); - EXPECT_DOUBLE_EQ(rate->activationEnergy_R(H_big_R), H_big_R); - EXPECT_DOUBLE_EQ(rate->activationEnergy_R(H_small_R), 0); - EXPECT_NEAR(rate->activationEnergy_R(H_mid_R), Ea / GasConstant, 1e-7); + EXPECT_DOUBLE_EQ(rate->effectiveActivationEnergy_R(H_big_R), H_big_R); + EXPECT_DOUBLE_EQ(rate->effectiveActivationEnergy_R(H_small_R), 0); + EXPECT_NEAR(rate->effectiveActivationEnergy_R(H_mid_R), Ea / GasConstant, 1e-7); EXPECT_TRUE(rate->allowNegativePreExponentialFactor()); EXPECT_FALSE(R->allow_negative_orders); } From 3441bc208d4d38ce9a083bceb1fcc6fc7f63da70 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 6 Feb 2022 16:29:22 -0600 Subject: [PATCH 06/10] [Kinetics] Define TwoTempPlasmaRate as BulkRate --- include/cantera/kinetics/Arrhenius.h | 64 ++++++++++------------------ src/kinetics/Arrhenius.cpp | 29 +++---------- 2 files changed, 29 insertions(+), 64 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 9eca8673d6..75e5f78cbb 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -210,12 +210,13 @@ class Arrhenius3 : public ArrheniusBase * Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures. * Plasma Sources Science and Technology, 1(3), 207. * doi: 10.1088/0963-0252/1/3/011 + * * @ingroup arrheniusGroup */ -class TwoTempPlasmaRate final : public ArrheniusBase, public ReactionRate +class TwoTempPlasma : public ArrheniusBase { public: - TwoTempPlasmaRate(); + TwoTempPlasma(); //! Constructor. /*! @@ -225,43 +226,19 @@ class TwoTempPlasmaRate final : public ArrheniusBase, public ReactionRate * @param Ea Activation energy in energy units [J/kmol] * @param EE Activation electron energy in energy units [J/kmol] */ - TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); - - unique_ptr newMultiRate() const override { - return unique_ptr( - new MultiRate); - } + TwoTempPlasma(double A, double b, double Ea=0.0, double EE=0.0); - //! Constructor based on AnyMap content - TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}) - : TwoTempPlasmaRate() - { - setParameters(node, rate_units); - } - - //! Identifier of reaction rate type - virtual const std::string type() const override { + const std::string rateType() const { return "two-temperature-plasma"; } - //! Perform object setup based on AnyMap node information - /*! - * @param node AnyMap containing rate information - * @param rate_units Unit definitions specific to rate information - */ - virtual void setParameters(const AnyMap& node, const UnitStack& rate_units) override; - - virtual void getParameters(AnyMap& node) const override; - - void check(const std::string& equation, const AnyMap& node) override { - checkRate(equation, node); + //! Context (unused) + void setRateContext(const Reaction& rxn, const Kinetics& kin) { } //! Evaluate reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double evalFromStruct(const TwoTempPlasmaData& shared_data) { + //! @internal Non-virtual method that should not be overloaded + double evalRate(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 + @@ -269,11 +246,14 @@ class TwoTempPlasmaRate final : public ArrheniusBase, public ReactionRate * shared_data.recipTe * shared_data.recipT); } - //! Return the activation energy divided by the gas constant (i.e. the - //! activation temperature) [K] - double activationEnergy_R() const { - return m_Ea_R; - } + //! Evaluate derivative of reaction rate with respect to temperature + //! divided by reaction rate + /*! + * This method does not consider changes of electron temperature. + * A corresponding warning is raised. + * @internal Non-virtual method that should not be overloaded + */ + double ddTScaled(const TwoTempPlasmaData& shared_data) const; //! Return the electron activation energy *Ea* [J/kmol] double activationElectronEnergy() const { @@ -338,6 +318,7 @@ class BlowersMasel : public ArrheniusBase return "Blowers-Masel"; } + //! Set context void setRateContext(const Reaction& rxn, const Kinetics& kin); //! Update information specific to reaction @@ -456,12 +437,12 @@ class BulkRate final : public RateType, public ReactionRate } void check(const std::string& equation, const AnyMap& node) override { - checkRate(equation, node); + RateType::checkRate(equation, node); } virtual void setContext(const Reaction& rxn, const Kinetics& kin) override { // as this method is virtual, it cannot be templated - setRateContext(rxn, kin); + RateType::setRateContext(rxn, kin); } //! Evaluate reaction rate @@ -469,7 +450,7 @@ class BulkRate final : public RateType, public ReactionRate * @param shared_data data shared by all reactions of a given type */ double evalFromStruct(const DataType& shared_data) const { - return evalRate(shared_data); + return RateType::evalRate(shared_data); } //! Evaluate derivative of reaction rate with respect to temperature @@ -478,11 +459,12 @@ class BulkRate final : public RateType, public ReactionRate * @param shared_data data shared by all reactions of a given type */ double ddTScaledFromStruct(const DataType& shared_data) const { - return ddTScaled(shared_data); + return RateType::ddTScaled(shared_data); } }; typedef BulkRate ArrheniusRate; +typedef BulkRate TwoTempPlasmaRate; typedef BulkRate BlowersMaselRate; } diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index c6319b6174..8aab83b4e7 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -126,14 +126,14 @@ void ArrheniusBase::checkRate(const std::string& equation, const AnyMap& node) } } -TwoTempPlasmaRate::TwoTempPlasmaRate() +TwoTempPlasma::TwoTempPlasma() : ArrheniusBase() { m_Ea_str = "Ea-gas"; m_E4_str = "Ea-electron"; } -TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) +TwoTempPlasma::TwoTempPlasma(double A, double b, double Ea, double EE) : ArrheniusBase(A, b, Ea) { m_Ea_str = "Ea-gas"; @@ -141,28 +141,11 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) m_E4_R = EE / GasConstant; } -void TwoTempPlasmaRate::setParameters(const AnyMap& node, const UnitStack& rate_units) +double TwoTempPlasma::ddTScaled(const TwoTempPlasmaData& shared_data) const { - m_negativeA_ok = node.getBool("negative-A", false); - if (!node.hasKey("rate-constant")) { - setRateParameters(AnyValue(), node.units(), rate_units); - return; - } - setRateParameters(node["rate-constant"], node.units(), rate_units); -} - -void TwoTempPlasmaRate::getParameters(AnyMap& rateNode) const -{ - if (m_negativeA_ok) { - rateNode["negative-A"] = true; - } - AnyMap node; - ArrheniusBase::getRateParameters(node); - if (!node.empty()) { - // object is configured - rateNode["rate-constant"] = std::move(node); - } - rateNode["type"] = type(); + warn_user("TwoTempPlasma::ddTScaled", + "Temperature derivative does not consider changes of electron temperature."); + return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; } BlowersMasel::BlowersMasel() From 62979a015caa9cccbe8e9baa8fe42531b85be3f4 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 7 Feb 2022 09:43:08 -0600 Subject: [PATCH 07/10] Eliminate C++ TwoTempPlasmaReaction class specialization --- include/cantera/kinetics/Arrhenius.h | 5 +-- include/cantera/kinetics/Reaction.h | 18 -------- interfaces/cython/cantera/_cantera.pxd | 3 -- interfaces/cython/cantera/reaction.pyx | 59 ++------------------------ src/kinetics/Arrhenius.cpp | 11 +++++ src/kinetics/Reaction.cpp | 35 --------------- src/kinetics/ReactionFactory.cpp | 5 +-- 7 files changed, 17 insertions(+), 119 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 75e5f78cbb..598a3212a1 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -232,9 +232,8 @@ class TwoTempPlasma : public ArrheniusBase return "two-temperature-plasma"; } - //! Context (unused) - void setRateContext(const Reaction& rxn, const Kinetics& kin) { - } + //! Context + void setRateContext(const Reaction& rxn, const Kinetics& kin); //! Evaluate reaction rate //! @internal Non-virtual method that should not be overloaded diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index e959785a0b..e759c13f39 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -528,24 +528,6 @@ class ThreeBodyReaction3 : public Reaction }; -//! A reaction for two-temperature-plasma reaction rate -class TwoTempPlasmaReaction: public Reaction -{ -public: - TwoTempPlasmaReaction(); - TwoTempPlasmaReaction(const Composition& reactants, const Composition& products, - const TwoTempPlasmaRate& rate); - - TwoTempPlasmaReaction(const AnyMap& node, const Kinetics& kin); - - virtual std::string type() const { - return "two-temperature-plasma"; - } - - virtual void validate(); -}; - - //! A falloff reaction that is first-order in [M] at low pressure, like a third-body //! reaction, but zeroth-order in [M] as pressure increases. //! In addition, the class supports chemically-activated reactions where the rate diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 1895268fb7..d8eebe07c5 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -666,9 +666,6 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxThreeBodyReaction3 "Cantera::ThreeBodyReaction3" (CxxReaction): CxxThreeBodyReaction3() - cdef cppclass CxxTwoTempPlasmaReaction "Cantera::TwoTempPlasmaReaction"(CxxReaction): - CxxTwoTempPlasmaReaction() - cdef cppclass CxxFalloffReaction3 "Cantera::FalloffReaction3" (CxxReaction): CxxFalloffReaction3() diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index 22190d73e6..ef8046c5cb 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -290,6 +290,9 @@ cdef class TwoTempPlasmaRate(ArrheniusTypeRate): def __cinit__(self, A=None, b=None, Ea_gas=0.0, Ea_electron=0.0, input_data=None, init=True): if init: + if A is None and b is None: + Ea_gas = None + Ea_electron = None self._cinit(input_data, A=A, b=b, Ea_gas=Ea_gas, Ea_electron=Ea_electron) def __call__(self, double temperature, double elec_temp): @@ -2414,62 +2417,6 @@ cdef wrapBlowersMasel(CxxBlowersMasel2* rate, Reaction reaction): return r -cdef class TwoTempPlasmaReaction(Reaction): - """ - A reaction which rate coefficient depends on both gas and electron temperature - - An example for the definition of an `TwoTempPlasmaReaction` object is given as:: - - rxn = TwoTempPlasmaReaction( - equation="O2 + E <=> O2-", - rate={"A": 17283, "b": -3.1, "Ea-gas": -5820088, "Ea-electron": 10808733}, - kinetics=gas) - - The YAML description corresponding to this reaction is:: - - equation: O2 + E <=> O2- - type: two-temperature-plasma - rate-constant: {A: 17283, b: -3.1, Ea-gas: -700 K, Ea-electron: 1300 K} - """ - _reaction_type = "two-temperature-plasma" - - def __init__(self, reactants=None, products=None, rate=None, *, equation=None, - Kinetics kinetics=None, init=True, **kwargs): - - if reactants and products and not equation: - equation = self.equation - - if init and equation and kinetics: - rxn_type = self._reaction_type - spec = {"equation": equation, "type": rxn_type} - if isinstance(rate, dict): - replaced_rate = {} - for key, value in rate.items(): - replaced_rate[key.replace("_", "-")] = value - spec["rate-constant"] = replaced_rate - elif rate is None or isinstance(rate, TwoTempPlasmaRate): - pass - else: - raise TypeError("Invalid rate definition") - - self._reaction = CxxNewReaction(dict_to_anymap(spec), - deref(kinetics.kinetics)) - self.reaction = self._reaction.get() - - if isinstance(rate, TwoTempPlasmaRate): - self.rate = rate - - property rate: - """ Get/Set the `TwoTempPlasmaRate` rate coefficients for this reaction. """ - def __get__(self): - cdef CxxTwoTempPlasmaReaction* r = self.reaction - return TwoTempPlasmaRate.wrap(r.rate()) - def __set__(self, rate): - cdef CxxTwoTempPlasmaReaction* r = self.reaction - cdef TwoTempPlasmaRate rate_ = rate - r.setRate(rate_._rate) - - cdef class InterfaceReaction(ElementaryReaction): """ A reaction occurring on an `Interface` (i.e. a surface or an edge) """ _reaction_type = "interface" diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index 8aab83b4e7..372961e07e 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -148,6 +148,17 @@ double TwoTempPlasma::ddTScaled(const TwoTempPlasmaData& shared_data) const return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; } +void TwoTempPlasma::setRateContext(const Reaction& rxn, const Kinetics& kin) +{ + // TwoTempPlasmaReaction is for a non-equilirium plasma, and the reverse rate + // cannot be calculated from the conventional thermochemistry. + // @todo implement the reversible rate for non-equilibrium plasma + if (rxn.reversible) { + throw InputFileError("TwoTempPlasma::setRateContext", rxn.input, + "TwoTempPlasmaRate does not support reversible reactions"); + } +} + BlowersMasel::BlowersMasel() : m_deltaH_R(0.) { diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index ee62390ea9..c382a84173 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -1138,41 +1138,6 @@ std::string ThreeBodyReaction3::productString() const } } -TwoTempPlasmaReaction::TwoTempPlasmaReaction() -{ - setRate(newReactionRate(type())); -} - -TwoTempPlasmaReaction::TwoTempPlasmaReaction( - const Composition& reactants, const Composition& products, - const TwoTempPlasmaRate& rate) - : Reaction(reactants, products) -{ - m_rate.reset(new TwoTempPlasmaRate(rate)); -} - -TwoTempPlasmaReaction::TwoTempPlasmaReaction(const AnyMap& node, const Kinetics& kin) -{ - if (!node.empty()) { - setParameters(node, kin); - setRate(newReactionRate(node, calculateRateCoeffUnits3(kin))); - } else { - setRate(newReactionRate(type())); - } -} - -void TwoTempPlasmaReaction::validate() -{ - Reaction::validate(); - // TwoTempPlasmaReaction is for a non-equilirium plasma, and the reverse rate - // cannot be calculated from the conventional thermochemistry. - // @todo implement the reversible rate for non-equilibrium plasma - if (reversible) { - throw InputFileError("Reaction::validate", input, - "TwoTempPlasmaReaction may only be given for irreversible reactions"); - } -} - FalloffReaction3::FalloffReaction3() : Reaction() { diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index ea4095db41..01c7fa36a4 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -131,10 +131,7 @@ ReactionFactory::ReactionFactory() return R; }); - // register electron-temprature reactions - reg("two-temperature-plasma", [](const AnyMap& node, const Kinetics& kin) { - return new TwoTempPlasmaReaction(node, kin); - }); + addAlias("reaction", "two-temperature-plasma"); addAlias("reaction", "Blowers-Masel"); From d610a2c5ff90f744a904bf18e21b9dcf91a612f1 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 7 Feb 2022 07:57:12 -0600 Subject: [PATCH 08/10] [UnitTests] Update tests and introduce pytest features --- .../cython/cantera/test/test_reaction.py | 198 +++++++----------- 1 file changed, 81 insertions(+), 117 deletions(-) diff --git a/interfaces/cython/cantera/test/test_reaction.py b/interfaces/cython/cantera/test/test_reaction.py index 5c5e4f19ab..138eb7e026 100644 --- a/interfaces/cython/cantera/test/test_reaction.py +++ b/interfaces/cython/cantera/test/test_reaction.py @@ -5,6 +5,7 @@ import cantera as ct import numpy as np from . import utilities +import pytest class TestImplicitThirdBody(utilities.CanteraTest): @@ -71,7 +72,7 @@ def test_duplicate(self): gas1.write_yaml(fname) with self.assertRaisesRegex(Exception, "Undeclared duplicate reactions"): - gas2 = ct.Solution(fname) + ct.Solution(fname) Path(fname).unlink() @@ -166,7 +167,7 @@ def from_yaml(self): def from_dict(self, input=None): # create reaction rate object from dictionary if input is None: - input = self.from_parts().input_data + input = self.from_yaml().input_data else: self.assertIsInstance(input, dict) return ct.ReactionRate.from_dict(input) @@ -210,7 +211,7 @@ def test_unconfigured(self): def test_roundtrip(self): # check round-trip instantiation via input_data - rate0 = self.from_parts() + rate0 = self.from_yaml() input_data = rate0.input_data rate1 = self.from_dict(input_data) self.assertNear(self.eval(rate0), self.eval(rate1)) @@ -222,7 +223,7 @@ def test_with_units(self): with self.assertRaisesRegex(Exception, "not supported"): ct.ReactionRate.from_yaml(yaml) - def test_temperature_derivative(self): + def test_derivative_ddT(self): # check temperature derivative against numerical derivative deltaT = self.gas.derivative_settings["rtol-delta"] deltaT *= self.gas.T @@ -237,7 +238,7 @@ def test_temperature_derivative(self): k1 = self.eval(rate) self.assertNear((k1 - k0) / deltaT, drate[self._index], 1e-6) - def test_pressure_derivative(self): + def test_derivative_ddP(self): # check pressure derivative against numerical derivative deltaP = self.gas.derivative_settings["rtol-delta"] deltaP *= self.gas.P @@ -257,13 +258,9 @@ class TestArrheniusRate(ReactionRateTests, utilities.CanteraTest): _type = "Arrhenius" _index = 0 _input = {"rate-constant": {"A": 38.7, "b": 2.7, "Ea": 26191840.0}} + _parts = {"A": 38.7, "b": 2.7, "Ea": 26191840.0} _yaml = "rate-constant: {A: 38.7, b: 2.7, Ea: 6260.0 cal/mol}" - @classmethod - def setUpClass(cls): - ReactionRateTests.setUpClass() - cls._parts = cls._input["rate-constant"] - def test_from_parts(self): rate = self.from_parts() self.assertEqual(self._parts["A"], rate.pre_exponential_factor) @@ -284,7 +281,7 @@ def test_standalone(self): with self.assertRaisesRegex(Exception, "not supported"): ct.ReactionRate.from_yaml(yaml) - def test_temperature_derivative_exact(self): + def test_derivative_ddT_exact(self): # check exact derivative against analytical and numerical derivatives rate = self.from_parts() T = 1000. @@ -314,16 +311,12 @@ class TestBlowersMaselRate(ReactionRateTests, utilities.CanteraTest): _type = "Blowers-Masel" _index = 6 _input = {"rate-constant": {"A": 38700, "b": 2.7, "Ea0": 1.0958665856e8, "w": 1.7505856e13}} + _parts = {"A": 38700, "b": 2.7, "Ea0": 1.0958665856e8, "w": 1.7505856e13} _yaml = """ type: Blowers-Masel rate-constant: {A: 38700, b: 2.7, Ea0: 1.0958665856e8, w: 1.7505856e13} """ - @classmethod - def setUpClass(cls): - ReactionRateTests.setUpClass() - cls._parts = cls._input["rate-constant"] - def eval(self, rate): delta_enthalpy = self.gas.delta_enthalpy[self._index] return rate(self.gas.T, delta_enthalpy) @@ -343,9 +336,9 @@ def test_negative_A(self): rate.allow_negative_pre_exponential_factor = True self.assertTrue(rate.allow_negative_pre_exponential_factor) - @utilities.unittest.skip("change of reaction enthalpy is not considered") - def test_temperature_derivative(self): - super().test_temperature_derivative() + @pytest.mark.skip("Change of reaction enthalpy is not considered") + def test_derivative_ddT(self): + super().test_derivative_ddT() class TestTwoTempPlasmaRate(ReactionRateTests, utilities.CanteraTest): @@ -355,16 +348,12 @@ class TestTwoTempPlasmaRate(ReactionRateTests, utilities.CanteraTest): _type = "two-temperature-plasma" _index = 11 _input = {"rate-constant": {"A": 17283, "b": -3.1, "Ea-gas": -5820000, "Ea-electron": 1081000}} + _parts = {"A": 17283, "b": -3.1, "Ea_gas": -5820000, "Ea_electron": 1081000} _yaml = """ type: two-temperature-plasma rate-constant: {A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol} """ - @classmethod - def setUpClass(cls): - ReactionRateTests.setUpClass() - cls._parts = {key.replace("-", "_"): value for key, value in cls._input["rate-constant"].items()} - def eval(self, rate): # check evaluation as a function of temperature and electron temperature return rate(self.gas.T, self.gas.Te) @@ -391,17 +380,13 @@ class TestTwoTempPlasmaRateShort(ReactionRateTests, utilities.CanteraTest): _cls = ct.TwoTempPlasmaRate _type = "two-temperature-plasma" _index = 12 - _input = {"rate-constant": {"A": 17283, "b": -3.1, "Ea-gas": 0.0, "Ea-electron": 0.0}} + _input = {"rate-constant": {"A": 17283, "b": -3.1}} + _parts = {"A": 17283, "b": -3.1} _yaml = """ type: two-temperature-plasma rate-constant: {A: 17283, b: -3.1, Ea-gas: 0.0 J/mol, Ea-electron: 0.0 J/mol} """ - @classmethod - def setUpClass(cls): - ReactionRateTests.setUpClass() - cls._parts = {key.replace("-", "_"): value for key, value in cls._input["rate-constant"].items()} - def eval(self, rate): # check evaluation as a function of temperature and electron temperature return rate(self.gas.T, self.gas.Te) @@ -410,17 +395,10 @@ def test_from_parts(self): rate = self.from_parts() self.assertEqual(self._parts["A"], rate.pre_exponential_factor) self.assertEqual(self._parts["b"], rate.temperature_exponent) - self.assertAlmostEqual(self._parts["Ea_gas"], rate.activation_energy) - self.assertAlmostEqual(self._parts["Ea_electron"], rate.activation_electron_energy) + self.assertAlmostEqual(rate.activation_energy, 0.) + self.assertAlmostEqual(rate.activation_electron_energy, 0.) self.check_rate(rate) - def test_negative_A(self): - # test reaction rate property - rate = self.from_parts() - self.assertFalse(rate.allow_negative_pre_exponential_factor) - rate.allow_negative_pre_exponential_factor = True - self.assertTrue(rate.allow_negative_pre_exponential_factor) - class FalloffRateTests(ReactionRateTests): # test Falloff rate expressions @@ -443,7 +421,7 @@ def test_data(self): for n in self._n_data: rate.falloff_coeffs = np.random.rand(n) - def test_temperature_derivative(self): + def test_derivative_ddT(self): pert = self.gas.derivative_settings["rtol-delta"] deltaT = self.gas.T * pert TP = self.gas.TP @@ -464,7 +442,7 @@ def test_temperature_derivative(self): k1 = self.eval(rate) self.assertNear((k1 - k0) / deltaT, drate[self._index], 1e-6) - def test_pressure_derivative(self): + def test_derivative_ddP(self): pert = self.gas.derivative_settings["rtol-delta"] deltaP = self.gas.P * pert rate = self.from_yaml() @@ -709,6 +687,7 @@ class ReactionTests: _rxn_type = None # name of reaction type _rate_type = None # name of reaction rate type _legacy = False # object uses legacy framework + _legacy_uses_rate = True # legacy object implements rate property _yaml = None # YAML parameterization _deprecated_getters = {} # test legacy getters (old framework) _deprecated_setters = {} # test legacy setters (old framework) @@ -740,21 +719,35 @@ def from_dict(self): input_data = self.from_yaml().input_data return ct.Reaction.from_dict(input_data, kinetics=self.gas) - def from_rate(self, rate): - if rate is None and self._rate_cls is not None: + def from_empty(self): + # create reaction object with uninitialized rate + if self._rate_cls is None: + rate = None + else: # Create an "empty" rate of the correct type for merged reaction types where # the only way they are distinguished is by the rate type - rate = self._rate_cls() - # create reaction object from keywords / rate + rate=self._rate_cls() + return self._cls(equation=self._equation, rate=rate, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) + + def from_rate(self, rate): + if not self._legacy_uses_rate: + pytest.skip("Legacy: rate object is not defined [1]") + if rate is None and self._rate is None: + # this does not work when no specialized reaction class exists + pytest.skip("Creation from dictionary is not supported") return self._cls(equation=self._equation, rate=rate, kinetics=self.gas, legacy=self._legacy, **self._kwargs) def from_parts(self): # create reaction rate object from parts + if self._rate_obj is None: + pytest.skip("Legacy: rate object is not defined [2]") orig = self.gas.reaction(self._index) if self._legacy: rxn = self._cls(orig.reactants, orig.products, legacy=self._legacy) rxn.rate = self._rate_obj + rxn.reversible = "<=>" in self._equation return rxn else: return self._cls(orig.reactants, orig.products, rate=self._rate_obj, @@ -799,41 +792,38 @@ def check_solution(self, gas2, check_legacy=True): def test_rate(self): # check consistency of reaction rate and forward rate constant - if self._rate_obj is not None: - self.check_rate(self._rate_obj) + if self._rate_obj is None: + pytest.skip("Legacy: rate object is not defined [3]") + self.check_rate(self._rate_obj) def test_from_rate(self): # check instantiation from keywords / rate defined by dictionary - if self._rate is not None: - self.check_rxn(self.from_rate(self._rate)) + self.check_rxn(self.from_rate(self._rate)) def test_from_rate_obj(self): # check instantiation from keywords / rate provided as object - if self._rate_obj is not None: - self.check_rxn(self.from_rate(self._rate_obj)) + self.check_rxn(self.from_rate(self._rate_obj)) def test_from_parts(self): # check instantiation from parts (reactants, products, rate expression) - if self._rate_obj is not None: - self.check_rxn(self.from_parts()) + self.check_rxn(self.from_parts()) def test_from_yaml(self): # check constructors (from yaml input) - if self._yaml is not None: - # check instantiation from yaml string - self.check_rxn(self.from_yaml()) - self.check_rxn(self.from_yaml(deprecated=True)) + self.check_rxn(self.from_yaml()) + self.check_rxn(self.from_yaml(deprecated=True)) def test_from_dict(self): # check instantiation from a yaml dictionary (input_data) - if self._yaml is not None: - # cannot compare types as input_data does not recreate legacy objects - self.check_rxn(self.from_dict(), check_legacy=False) + # cannot compare types as input_data does not recreate legacy objects + self.check_rxn(self.from_dict(), check_legacy=False) def test_add_rxn(self): # check adding new reaction to solution - if self._rate_obj is None: - return + if self._yaml is not None: + rxn = self.from_yaml() + else: + rxn = self.from_rate(self._rate_obj) gas2 = ct.Solution(thermo="IdealGas", kinetics="GasKinetics", species=self.species, reactions=[]) gas2.TPX = self.gas.TPX @@ -856,7 +846,7 @@ def test_no_rate(self): # check behavior for instantiation from keywords / no rate if self._rate_obj is None: return - rxn = self.from_rate(None) + rxn = self.from_empty() if self._legacy: self.assertNear(rxn.rate(self.gas.T), 0.) else: @@ -877,16 +867,21 @@ def test_no_rate(self): def test_replace_rate(self): # check replacing reaction rate expression - if self._rate_obj is None: - return - rxn = self.from_rate(None) + if not self._legacy_uses_rate: + pytest.skip("Legacy: rate property not implemented") + elif self._yaml is not None: + rxn = self.from_yaml() + else: + rxn = self.from_rate(self._rate_obj) + rxn = self.from_empty() rxn.rate = self._rate_obj self.check_rxn(rxn) def test_roundtrip(self): # check round-trip instantiation via input_data + rxn = self.from_yaml() if self._legacy: - return + pytest.skip("Legacy: round-trip conversion is not supported") rxn = self.from_rate(self._rate_obj) rate_input_data = rxn.rate.input_data rate_obj = rxn.rate.__class__(input_data=rate_input_data) @@ -903,9 +898,6 @@ def check_equal(self, one, two): def test_deprecated_getters(self): # check property getters deprecated in new framework - if self._yaml is None: - return - rxn = self.from_yaml() for attr, default in self._deprecated_getters.items(): if self._legacy: @@ -920,9 +912,6 @@ def test_deprecated_getters(self): def test_deprecated_setters(self): # check property setters deprecated in new framework - if self._yaml is None: - return - rxn = self.from_yaml() for attr, new in self._deprecated_setters.items(): if self._legacy: @@ -936,9 +925,6 @@ def test_deprecated_setters(self): def test_deprecated_callers(self): # check methods deprecated in new framework - if self._yaml is None: - return - rxn = self.from_yaml() for state, value in self._deprecated_callers.items(): T, P = state @@ -978,7 +964,7 @@ def setUpClass(cls): def test_arrhenius(self): # test assigning Arrhenius rate rate = ct.Arrhenius(self._rate["A"], self._rate["b"], self._rate["Ea"]) - rxn = self.from_rate(None) + rxn = self.from_empty() if self._legacy: rxn.rate = rate else: @@ -1023,10 +1009,6 @@ def from_parts(self): rxn.efficiencies = self._kwargs["efficiencies"] return rxn - def test_rate(self): - # rate constant contains third-body concentration - pass - def test_efficiencies(self): # check efficiencies rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas, @@ -1077,11 +1059,12 @@ def test_efficiencies(self): class TestTwoTempPlasma(ReactionTests, utilities.CanteraTest): # test two-temperature plasma reaction - _cls = ct.TwoTempPlasmaReaction - _rxn_type = "two-temperature-plasma" + _cls = ct.Reaction + _rate_cls = ct.TwoTempPlasmaRate + _rxn_type = "reaction" _rate_type = "two-temperature-plasma" _equation = "O + H => O + H" - _rate = {"A": 17283, "b": -3.1, "Ea_gas": -5820000, "Ea_electron": 1081000} + _rate_obj = ct.TwoTempPlasmaRate(A=17283, b=-3.1, Ea_gas=-5820000, Ea_electron=1081000) _index = 11 _yaml = """ equation: O + H => O + H @@ -1089,22 +1072,19 @@ class TestTwoTempPlasma(ReactionTests, utilities.CanteraTest): rate-constant: {A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol} """ - @classmethod - def setUpClass(cls): - ReactionTests.setUpClass() - cls._rate_obj = ct.TwoTempPlasmaRate(**cls._rate) - def from_parts(self): - # create reaction rate object from parts - orig = self.gas.reaction(self._index) - rxn = self._cls(orig.reactants, orig.products, legacy=self._legacy) - rxn.rate = self._rate_obj + rxn = super().from_parts() rxn.reversible = False return rxn def eval_rate(self, rate): return rate(self.gas.T, self.gas.Te) + def test_reversible(self): + rxn = super().from_parts() + with self.assertRaisesRegex(ct.CanteraError, "does not support reversible"): + self.check_rxn(rxn) + class TestBlowersMasel(ReactionTests, utilities.CanteraTest): # test updated version of Blowers-Masel reaction @@ -1136,6 +1116,7 @@ class TestTroe2(ReactionTests, utilities.CanteraTest): _index = 2 _rxn_type = "falloff-legacy" _legacy = True + _legacy_uses_rate = False _yaml = """ equation: 2 OH (+ M) <=> H2O2 (+ M) # Reaction 3 type: falloff-legacy @@ -1145,19 +1126,6 @@ class TestTroe2(ReactionTests, utilities.CanteraTest): efficiencies: {AR: 0.7, H2: 2.0, H2O: 6.0} """ - def from_parts(self): - rxn = ReactionTests.from_parts(self) - rxn.efficiencies = self._kwargs["efficiencies"] - return rxn - - def test_from_rate(self): - # do not port creation from legacy Fallout objects - pass - - def test_from_yaml(self): - # check constructors (from yaml input) - self.check_rxn(self.from_yaml(), check_legacy=False) - class TestTroe(ReactionTests, utilities.CanteraTest): @@ -1254,6 +1222,7 @@ class TestChemicallyActivated2(ReactionTests, utilities.CanteraTest): _index = 10 _rxn_type = "chemically-activated-legacy" _legacy = True + _legacy_uses_rate = False _yaml = """ equation: H2O + OH (+M) <=> HO2 + H2 (+M) # Reaction 11 units: {length: cm, quantity: mol, activation-energy: cal/mol} @@ -1262,14 +1231,6 @@ class TestChemicallyActivated2(ReactionTests, utilities.CanteraTest): high-P-rate-constant: [5.88E-14, 6.721, -3022.227] """ - def test_from_rate(self): - # do not port creation from legacy Fallout objects - pass - - def test_from_yaml(self): - # check constructors (from yaml input) - self.check_rxn(self.from_yaml(), check_legacy=False) - class TestChemicallyActivated(ReactionTests, utilities.CanteraTest): # test Chemically Activated falloff reaction @@ -1319,6 +1280,7 @@ class TestPlog2(ReactionTests, utilities.CanteraTest): _index = 3 _rxn_type = "pressure-dependent-Arrhenius-legacy" _legacy = True + _legacy_uses_rate = False _yaml = """ equation: H2 + O2 <=> 2 OH type: pressure-dependent-Arrhenius-legacy @@ -1380,6 +1342,7 @@ class TestPlog(TestPlog2): _rxn_type = "reaction" _rate_type = "pressure-dependent-Arrhenius" _legacy = False + _legacy_uses_rate = True _yaml = """ equation: H2 + O2 <=> 2 OH type: pressure-dependent-Arrhenius @@ -1411,6 +1374,7 @@ class TestChebyshev2(ReactionTests, utilities.CanteraTest): _index = 4 _rxn_type = "Chebyshev-legacy" _legacy = True + _legacy_uses_rate = False _yaml = """ equation: HO2 <=> OH + O type: Chebyshev-legacy @@ -1447,6 +1411,7 @@ class TestChebyshev(TestChebyshev2): [ 1.9764e+00, 1.0037e+00, 7.2865e-03, -3.0432e-02], [ 3.1770e-01, 2.6889e-01, 9.4806e-02, -7.6385e-03]]) _legacy = False + _legacy_uses_rate = True _yaml = """ equation: HO2 <=> OH + O type: Chebyshev @@ -1480,9 +1445,8 @@ def setUp(self): super().setUp() self._rate = lambda T: 38.7 * T**2.7 * exp(-3150.15428/T) - def test_roundtrip(self): - # overload default tester for round trip - pass + def from_yaml(self, deprecated=False): + pytest.skip("CustomReaction does not support YAML") def test_raises_invalid_rate(self): # check exception for instantiation from keywords / invalid rate From 49ff30b75f95e2dd80ed728477b0efd5c944d9a0 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 7 Feb 2022 12:46:21 -0600 Subject: [PATCH 09/10] [Kinetics] Eliminate multiple inheritance While multiple inheritance avoids the introduction of virtual functions, a direct inheritance scheme may avoid a small speed penalty. --- include/cantera/kinetics/Arrhenius.h | 38 +++++++------------------ include/cantera/kinetics/ReactionRate.h | 6 +++- include/cantera/kinetics/RxnRates.h | 6 +++- interfaces/cython/cantera/_cantera.pxd | 22 +++++++------- interfaces/cython/cantera/reaction.pyx | 9 ++---- src/kinetics/Arrhenius.cpp | 12 ++++---- src/kinetics/Falloff.cpp | 8 +++--- 7 files changed, 44 insertions(+), 57 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 598a3212a1..0b9e63e9a5 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -31,7 +31,7 @@ class AnyMap; * access from derived classes as well as classes that use Arrhenius-type expressions * internally, for example FalloffRate and PlogRate. */ -class ArrheniusBase +class ArrheniusBase : public ReactionRate { public: //! Default constructor. @@ -67,7 +67,7 @@ class ArrheniusBase void getRateParameters(AnyMap& node) const; //! Check rate expression - void checkRate(const std::string& equation, const AnyMap& node); + virtual void check(const std::string& equation, const AnyMap& node) override; //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending //! on the reaction order) @@ -161,14 +161,10 @@ class Arrhenius3 : public ArrheniusBase public: using ArrheniusBase::ArrheniusBase; // inherit constructors - const std::string rateType() const { + virtual const std::string type() const override { return "Arrhenius"; } - //! Context (unused) - void setRateContext(const Reaction& rxn, const Kinetics& kin) { - } - //! Evaluate reaction rate //! @internal Non-virtual method that should not be overloaded double evalRate(double logT, double recipT) const { @@ -228,12 +224,12 @@ class TwoTempPlasma : public ArrheniusBase */ TwoTempPlasma(double A, double b, double Ea=0.0, double EE=0.0); - const std::string rateType() const { + virtual const std::string type() const override { return "two-temperature-plasma"; } //! Context - void setRateContext(const Reaction& rxn, const Kinetics& kin); + virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; //! Evaluate reaction rate //! @internal Non-virtual method that should not be overloaded @@ -313,15 +309,14 @@ class BlowersMasel : public ArrheniusBase */ BlowersMasel(double A, double b, double Ea0, double w); - const std::string rateType() const { + virtual const std::string type() const override { return "Blowers-Masel"; } //! Set context - void setRateContext(const Reaction& rxn, const Kinetics& kin); + virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; //! Update information specific to reaction - //! @todo Move to BulkRate (as template function via std::enable_if) void updateFromStruct(const BlowersMaselData& shared_data) { if (shared_data.ready) { m_deltaH_R = 0.; @@ -389,7 +384,7 @@ class BlowersMasel : public ArrheniusBase //! A class template for bulk phase reaction rate specifications template -class BulkRate final : public RateType, public ReactionRate +class BulkRate final : public RateType { public: BulkRate() = default; @@ -405,10 +400,6 @@ class BulkRate final : public RateType, public ReactionRate new MultiRate, DataType>); } - virtual const std::string type() const override { - return RateType::rateType(); - } - virtual void setParameters( const AnyMap& node, const UnitStack& rate_units) override { @@ -430,20 +421,11 @@ class BulkRate final : public RateType, public ReactionRate // RateType object is configured node["rate-constant"] = std::move(rateNode); } - if (RateType::rateType() != "Arrhenius") { - node["type"] = type(); + if (RateType::type() != "Arrhenius") { + node["type"] = RateType::type(); } } - void check(const std::string& equation, const AnyMap& node) override { - RateType::checkRate(equation, node); - } - - virtual void setContext(const Reaction& rxn, const Kinetics& kin) override { - // as this method is virtual, it cannot be templated - RateType::setRateContext(rxn, kin); - } - //! Evaluate reaction rate /*! * @param shared_data data shared by all reactions of a given type diff --git a/include/cantera/kinetics/ReactionRate.h b/include/cantera/kinetics/ReactionRate.h index ca722e5780..9022768167 100644 --- a/include/cantera/kinetics/ReactionRate.h +++ b/include/cantera/kinetics/ReactionRate.h @@ -74,8 +74,12 @@ class ReactionRate //! //! where `RateType` is the derived class name and `DataType` is the corresponding //! container for parameters needed to evaluate reactions of that type. - virtual unique_ptr newMultiRate() const = 0; + virtual unique_ptr newMultiRate() const { + throw NotImplementedError("ReactionRate::newMultiRate", + "Not implemented by '{}' object.", type()); + } + //! String identifying reaction rate specialization virtual const std::string type() const = 0; //! Set parameters diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 4d530120b7..dce276a698 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -36,7 +36,7 @@ class Func1; * k_f = A T^b \exp (-E/RT) * \f] */ -class Arrhenius2 : public Arrhenius3 +class Arrhenius2 final : public Arrhenius3 { public: //! Default constructor. @@ -91,6 +91,10 @@ class Arrhenius2 : public Arrhenius3 doublereal updateRC(doublereal logT, doublereal recipT) const { return m_A * std::exp(m_b*logT - m_Ea_R*recipT); } + + virtual const std::string type() const override { + return "Arrhenius2"; + } }; diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index d8eebe07c5..890b127c47 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -457,7 +457,14 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef vector[shared_ptr[CxxReaction]] CxxGetReactions "getReactions" (XML_Node&) except +translate_exception cdef vector[shared_ptr[CxxReaction]] CxxGetReactions "getReactions" (CxxAnyValue&, CxxKinetics&) except +translate_exception - cdef cppclass CxxArrheniusBase "Cantera::ArrheniusBase": + cdef cppclass CxxReactionRate "Cantera::ReactionRate": + CxxReactionRate() + string type() + double eval(double) except +translate_exception + double eval(double, double) except +translate_exception + CxxAnyMap parameters() except +translate_exception + + cdef cppclass CxxArrheniusBase "Cantera::ArrheniusBase" (CxxReactionRate): CxxArrheniusBase() double preExponentialFactor() double temperatureExponent() @@ -472,23 +479,16 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxArrhenius2 "Cantera::Arrhenius2" (CxxArrhenius): CxxArrhenius2(double, double, double) - cdef cppclass CxxReactionRate "Cantera::ReactionRate": - CxxReactionRate() - string type() - double eval(double) except +translate_exception - double eval(double, double) except +translate_exception - CxxAnyMap parameters() except +translate_exception - - cdef cppclass CxxArrheniusRate "Cantera::ArrheniusRate" (CxxReactionRate, CxxArrhenius): + cdef cppclass CxxArrheniusRate "Cantera::ArrheniusRate" (CxxArrhenius): CxxArrheniusRate(CxxAnyMap) except +translate_exception CxxArrheniusRate(double, double, double) - cdef cppclass CxxTwoTempPlasmaRate "Cantera::TwoTempPlasmaRate" (CxxReactionRate, CxxArrheniusBase): + cdef cppclass CxxTwoTempPlasmaRate "Cantera::TwoTempPlasmaRate" (CxxArrheniusBase): CxxTwoTempPlasmaRate(CxxAnyMap) except +translate_exception CxxTwoTempPlasmaRate(double, double, double, double) double activationElectronEnergy() - cdef cppclass CxxBlowersMaselRate "Cantera::BlowersMaselRate" (CxxReactionRate, CxxArrheniusBase): + cdef cppclass CxxBlowersMaselRate "Cantera::BlowersMaselRate" (CxxArrheniusBase): CxxBlowersMaselRate(CxxAnyMap) except +translate_exception CxxBlowersMaselRate(double, double, double, double) double effectiveActivationEnergy(double) diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index ef8046c5cb..ae25faa160 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -214,8 +214,7 @@ cdef class ArrheniusRate(ArrheniusTypeRate): cdef set_cxx_object(self): self.rate = self._rate.get() - # CxxArrhenius does not have a common base with CxxReactionRate - self.base = self.cxx_object() + self.base = self.rate cdef CxxArrheniusRate* cxx_object(self): return self.rate @@ -248,8 +247,7 @@ cdef class BlowersMaselRate(ArrheniusTypeRate): cdef set_cxx_object(self): self.rate = self._rate.get() - # CxxArrheniusBase does not have a common base with CxxReactionRate - self.base = self.cxx_object() + self.base = self.rate cdef CxxBlowersMaselRate* cxx_object(self): return self.rate @@ -309,8 +307,7 @@ cdef class TwoTempPlasmaRate(ArrheniusTypeRate): cdef set_cxx_object(self): self.rate = self._rate.get() - # CxxArrheniusBase does not have a common base with CxxReactionRate - self.base = self.cxx_object() + self.base = self.rate cdef CxxTwoTempPlasmaRate* cxx_object(self): return self.rate diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index 372961e07e..4557c837a6 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -111,16 +111,16 @@ void ArrheniusBase::getRateParameters(AnyMap& node) const } -void ArrheniusBase::checkRate(const std::string& equation, const AnyMap& node) +void ArrheniusBase::check(const std::string& equation, const AnyMap& node) { if (!m_negativeA_ok && m_A < 0) { if (equation == "") { - throw CanteraError("ArrheniusBase::checkRate", + throw CanteraError("ArrheniusBase::check", "Detected negative pre-exponential factor (A={}).\n" "Enable 'allowNegativePreExponentialFactor' to suppress " "this message.", m_A); } - throw InputFileError("ArrheniusBase::checkRate", node, + throw InputFileError("ArrheniusBase::check", node, "Undeclared negative pre-exponential factor found in reaction '{}'", equation); } @@ -148,13 +148,13 @@ double TwoTempPlasma::ddTScaled(const TwoTempPlasmaData& shared_data) const return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; } -void TwoTempPlasma::setRateContext(const Reaction& rxn, const Kinetics& kin) +void TwoTempPlasma::setContext(const Reaction& rxn, const Kinetics& kin) { // TwoTempPlasmaReaction is for a non-equilirium plasma, and the reverse rate // cannot be calculated from the conventional thermochemistry. // @todo implement the reversible rate for non-equilibrium plasma if (rxn.reversible) { - throw InputFileError("TwoTempPlasma::setRateContext", rxn.input, + throw InputFileError("TwoTempPlasma::setContext", rxn.input, "TwoTempPlasmaRate does not support reversible reactions"); } } @@ -183,7 +183,7 @@ double BlowersMasel::ddTScaled(const BlowersMaselData& shared_data) const return (Ea_R * shared_data.recipT + m_b) * shared_data.recipT; } -void BlowersMasel::setRateContext(const Reaction& rxn, const Kinetics& kin) +void BlowersMasel::setContext(const Reaction& rxn, const Kinetics& kin) { m_stoich_coeffs.clear(); for (const auto& sp : rxn.reactants) { diff --git a/src/kinetics/Falloff.cpp b/src/kinetics/Falloff.cpp index ec5d18ba53..e06e5862b9 100644 --- a/src/kinetics/Falloff.cpp +++ b/src/kinetics/Falloff.cpp @@ -25,7 +25,7 @@ void FalloffRate::setLowRate(const Arrhenius3& low) { Arrhenius3 _low = low; _low.setAllowNegativePreExponentialFactor(m_negativeA_ok); - _low.checkRate("", AnyMap()); + _low.check("", AnyMap()); if (_low.preExponentialFactor() * m_highRate.preExponentialFactor() < 0.) { throw CanteraError("FalloffRate::setLowRate", "Detected inconsistent rate definitions;\nhigh and low " @@ -38,7 +38,7 @@ void FalloffRate::setHighRate(const Arrhenius3& high) { Arrhenius3 _high = high; _high.setAllowNegativePreExponentialFactor(m_negativeA_ok); - _high.checkRate("", AnyMap()); + _high.check("", AnyMap()); if (m_lowRate.preExponentialFactor() * _high.preExponentialFactor() < 0.) { throw CanteraError("FalloffRate::setHighRate", "Detected inconsistent rate definitions;\nhigh and low " @@ -118,8 +118,8 @@ void FalloffRate::getParameters(AnyMap& node) const void FalloffRate::check(const std::string& equation, const AnyMap& node) { - m_lowRate.checkRate(equation, node); - m_highRate.checkRate(equation, node); + m_lowRate.check(equation, node); + m_highRate.check(equation, node); double lowA = m_lowRate.preExponentialFactor(); double highA = m_highRate.preExponentialFactor(); From e0095a43a8c88581b28aee422d37cbe1c7a08553 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 7 Feb 2022 18:40:39 -0600 Subject: [PATCH 10/10] [Kinetics] Eliminate chained function calls --- include/cantera/kinetics/Arrhenius.h | 54 +++++++++++----------------- src/kinetics/Arrhenius.cpp | 8 ++--- 2 files changed, 25 insertions(+), 37 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 0b9e63e9a5..2ea4abaa1c 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -166,25 +166,29 @@ class Arrhenius3 : public ArrheniusBase } //! Evaluate reaction rate - //! @internal Non-virtual method that should not be overloaded double evalRate(double logT, double recipT) const { return m_A * std::exp(m_b * logT - m_Ea_R * recipT); } //! Evaluate natural logarithm of the rate constant. - //! @internal Non-virtual method that should not be overloaded double evalLog(double logT, double recipT) const { return m_logA + m_b * logT - m_Ea_R * recipT; } //! Evaluate reaction rate - double evalRate(const ArrheniusData& shared_data) const { + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const ArrheniusData& shared_data) const { return m_A * std::exp(m_b * shared_data.logT - m_Ea_R * shared_data.recipT); } //! Evaluate derivative of reaction rate with respect to temperature //! divided by reaction rate - double ddTScaled(const ArrheniusData& shared_data) const { + /*! + * @param shared_data data shared by all reactions of a given type + */ + double ddTScaledFromStruct(const ArrheniusData& shared_data) const { return (m_Ea_R * shared_data.recipT + m_b) * shared_data.recipT; } }; @@ -228,12 +232,13 @@ class TwoTempPlasma : public ArrheniusBase return "two-temperature-plasma"; } - //! Context virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; //! Evaluate reaction rate - //! @internal Non-virtual method that should not be overloaded - double evalRate(const TwoTempPlasmaData& shared_data) const { + /*! + * @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_b * shared_data.logTe - m_Ea_R * shared_data.recipT + @@ -246,9 +251,9 @@ class TwoTempPlasma : public ArrheniusBase /*! * This method does not consider changes of electron temperature. * A corresponding warning is raised. - * @internal Non-virtual method that should not be overloaded + * @param shared_data data shared by all reactions of a given type */ - double ddTScaled(const TwoTempPlasmaData& shared_data) const; + double ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const; //! Return the electron activation energy *Ea* [J/kmol] double activationElectronEnergy() const { @@ -313,7 +318,6 @@ class BlowersMasel : public ArrheniusBase return "Blowers-Masel"; } - //! Set context virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; //! Update information specific to reaction @@ -330,20 +334,21 @@ class BlowersMasel : public ArrheniusBase } //! Evaluate reaction rate - //! @internal Non-virtual method that should not be overloaded - double evalRate(const BlowersMaselData& shared_data) const { + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const BlowersMaselData& shared_data) const { double Ea_R = effectiveActivationEnergy_R(m_deltaH_R); return m_A * std::exp(m_b * shared_data.logT - Ea_R * shared_data.recipT); } - //! Evaluate derivative of reaction rate with respect to temperature //! divided by reaction rate /*! * This method does not consider potential changes due to a changed reaction * enthalpy. A corresponding warning is raised. - * @internal Non-virtual method that should not be overloaded + * @param shared_data data shared by all reactions of a given type */ - double ddTScaled(const BlowersMaselData& shared_data) const; + double ddTScaledFromStruct(const BlowersMaselData& shared_data) const; //! Return the effective activation energy (a function of the delta H of reaction) //! divided by the gas constant (i.e. the activation temperature) [K] @@ -384,7 +389,7 @@ class BlowersMasel : public ArrheniusBase //! A class template for bulk phase reaction rate specifications template -class BulkRate final : public RateType +class BulkRate : public RateType { public: BulkRate() = default; @@ -425,23 +430,6 @@ class BulkRate final : public RateType node["type"] = RateType::type(); } } - - //! Evaluate reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double evalFromStruct(const DataType& shared_data) const { - return RateType::evalRate(shared_data); - } - - //! Evaluate derivative of reaction rate with respect to temperature - //! divided by reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double ddTScaledFromStruct(const DataType& shared_data) const { - return RateType::ddTScaled(shared_data); - } }; typedef BulkRate ArrheniusRate; diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index 4557c837a6..469bba1171 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -141,16 +141,16 @@ TwoTempPlasma::TwoTempPlasma(double A, double b, double Ea, double EE) m_E4_R = EE / GasConstant; } -double TwoTempPlasma::ddTScaled(const TwoTempPlasmaData& shared_data) const +double TwoTempPlasma::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const { - warn_user("TwoTempPlasma::ddTScaled", + warn_user("TwoTempPlasma::ddTScaledFromStruct", "Temperature derivative does not consider changes of electron temperature."); return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; } void TwoTempPlasma::setContext(const Reaction& rxn, const Kinetics& kin) { - // TwoTempPlasmaReaction is for a non-equilirium plasma, and the reverse rate + // TwoTempPlasmaReaction is for a non-equilibrium plasma, and the reverse rate // cannot be calculated from the conventional thermochemistry. // @todo implement the reversible rate for non-equilibrium plasma if (rxn.reversible) { @@ -175,7 +175,7 @@ BlowersMasel::BlowersMasel(double A, double b, double Ea0, double w) m_E4_R = w / GasConstant; } -double BlowersMasel::ddTScaled(const BlowersMaselData& shared_data) const +double BlowersMasel::ddTScaledFromStruct(const BlowersMaselData& shared_data) const { warn_user("BlowersMasel::ddTScaledFromStruct", "Temperature derivative does not consider changes of reaction enthalpy.");