diff --git a/doc/sphinx/cython/kinetics.rst b/doc/sphinx/cython/kinetics.rst index fd88614160..c3ac82816c 100644 --- a/doc/sphinx/cython/kinetics.rst +++ b/doc/sphinx/cython/kinetics.rst @@ -149,6 +149,11 @@ StickingBlowersMaselRate Auxiliary Reaction Data ----------------------- +ThirdBody +^^^^^^^^^ +.. autoclass:: ThirdBody + :no-undoc-members: + Arrhenius ^^^^^^^^^ .. autoclass:: Arrhenius(A, b, E) diff --git a/include/cantera/kinetics/Falloff.h b/include/cantera/kinetics/Falloff.h index ff90d551d4..c982b4aba0 100644 --- a/include/cantera/kinetics/Falloff.h +++ b/include/cantera/kinetics/Falloff.h @@ -169,8 +169,11 @@ class FalloffRate : public ReactionRate return 0; } - virtual const std::string type() const { - return "Falloff"; + virtual const std::string type() const override { + if (m_chemicallyActivated) { + return "chemically-activated"; + } + return "falloff"; } //! Returns the number of parameters used by this parameterization. The @@ -179,7 +182,8 @@ class FalloffRate : public ReactionRate return 0; } - virtual void setParameters(const AnyMap& node, const UnitStack& rate_units); + virtual void setParameters( + const AnyMap& node, const UnitStack& rate_units) override; //! Get the values of the parameters for this object. *params* must be an //! array of at least nParameters() elements. @@ -220,7 +224,7 @@ class FalloffRate : public ReactionRate } virtual void check(const std::string& equation) override; - virtual void validate(const std::string& equation, const Kinetics& kin); + virtual void validate(const std::string& equation, const Kinetics& kin) override; //! Get flag indicating whether negative A values are permitted bool allowNegativePreExponentialFactor() const { @@ -297,12 +301,12 @@ class LindemannRate final : public FalloffRate setFalloffCoeffs(c); } - unique_ptr newMultiRate() const { + unique_ptr newMultiRate() const override{ return unique_ptr( new MultiRate); } - virtual const std::string type() const { + virtual const std::string subType() const override { return "Lindemann"; } }; @@ -358,7 +362,7 @@ class TroeRate final : public FalloffRate setFalloffCoeffs(c); } - unique_ptr newMultiRate() const { + unique_ptr newMultiRate() const override { return unique_ptr(new MultiRate); } @@ -367,9 +371,9 @@ class TroeRate final : public FalloffRate * @param c Vector of three or four doubles: The doubles are the parameters, * a, T_3, T_1, and (optionally) T_2 of the Troe parameterization */ - virtual void setFalloffCoeffs(const vector_fp& c); + virtual void setFalloffCoeffs(const vector_fp& c) override; - virtual void getFalloffCoeffs(vector_fp& c) const; + virtual void getFalloffCoeffs(vector_fp& c) const override; //! Update the temperature parameters in the representation /*! @@ -377,31 +381,32 @@ class TroeRate final : public FalloffRate * @param work Vector of working space, length 1, representing the * temperature-dependent part of the parameterization. */ - virtual void updateTemp(double T, double* work) const; + virtual void updateTemp(double T, double* work) const override; - virtual double F(double pr, const double* work) const; + virtual double F(double pr, const double* work) const override; - virtual size_t workSize() const { + virtual size_t workSize() const override { return 1; } - virtual const std::string type() const { + virtual const std::string subType() const override { return "Troe"; } - virtual size_t nParameters() const { + virtual size_t nParameters() const override { return 4; } - virtual void setParameters(const AnyMap& node, const UnitStack& rate_units); + virtual void setParameters( + const AnyMap& node, const UnitStack& rate_units) override; //! Sets params to contain, in order, \f[ (A, T_3, T_1, T_2) \f] /** * @deprecated To be removed after Cantera 3.0; superseded by getFalloffCoeffs() */ - virtual void getParameters(double* params) const; + virtual void getParameters(double* params) const override; - virtual void getParameters(AnyMap& node) const; + virtual void getParameters(AnyMap& node) const override; protected: //! parameter a in the 4-parameter Troe falloff function. Dimensionless @@ -460,7 +465,7 @@ class SriRate final : public FalloffRate setFalloffCoeffs(c); } - unique_ptr newMultiRate() const { + unique_ptr newMultiRate() const override { return unique_ptr(new MultiRate); } @@ -470,9 +475,9 @@ class SriRate final : public FalloffRate * a, b, c, d (optional; default 1.0), and e (optional; default * 0.0) of the SRI parameterization */ - virtual void setFalloffCoeffs(const vector_fp& c); + virtual void setFalloffCoeffs(const vector_fp& c) override; - virtual void getFalloffCoeffs(vector_fp& c) const; + virtual void getFalloffCoeffs(vector_fp& c) const override; //! Update the temperature parameters in the representation /*! @@ -480,31 +485,32 @@ class SriRate final : public FalloffRate * @param work Vector of working space, length 2, representing the * temperature-dependent part of the parameterization. */ - virtual void updateTemp(double T, double* work) const; + virtual void updateTemp(double T, double* work) const override; - virtual double F(double pr, const double* work) const; + virtual double F(double pr, const double* work) const override; - virtual size_t workSize() const { + virtual size_t workSize() const override { return 2; } - virtual const std::string type() const { + virtual const std::string subType() const override { return "SRI"; } - virtual size_t nParameters() const { + virtual size_t nParameters() const override { return 5; } - virtual void setParameters(const AnyMap& node, const UnitStack& rate_units); + virtual void setParameters( + const AnyMap& node, const UnitStack& rate_units) override; //! Sets params to contain, in order, \f[ (a, b, c, d, e) \f] /** * @deprecated To be removed after Cantera 3.0; superseded by getFalloffCoeffs() */ - virtual void getParameters(double* params) const; + virtual void getParameters(double* params) const override; - virtual void getParameters(AnyMap& node) const; + virtual void getParameters(AnyMap& node) const override; protected: //! parameter a in the 5-parameter SRI falloff function. Dimensionless. @@ -570,7 +576,7 @@ class TsangRate final : public FalloffRate setFalloffCoeffs(c); } - unique_ptr newMultiRate() const { + unique_ptr newMultiRate() const override { return unique_ptr(new MultiRate); } @@ -579,9 +585,9 @@ class TsangRate final : public FalloffRate * @param c Vector of one or two doubles: The doubles are the parameters, * a and (optionally) b of the Tsang F_cent parameterization */ - virtual void setFalloffCoeffs(const vector_fp& c); + virtual void setFalloffCoeffs(const vector_fp& c) override; - virtual void getFalloffCoeffs(vector_fp& c) const; + virtual void getFalloffCoeffs(vector_fp& c) const override; //! Update the temperature parameters in the representation /*! @@ -589,31 +595,32 @@ class TsangRate final : public FalloffRate * @param work Vector of working space, length 1, representing the * temperature-dependent part of the parameterization. */ - virtual void updateTemp(double T, double* work) const; + virtual void updateTemp(double T, double* work) const override; - virtual double F(double pr, const double* work) const; + virtual double F(double pr, const double* work) const override; - virtual size_t workSize() const { + virtual size_t workSize() const override { return 1; } - virtual const std::string type() const { + virtual const std::string subType() const override { return "Tsang"; } - virtual size_t nParameters() const { + virtual size_t nParameters() const override { return 2; } - virtual void setParameters(const AnyMap& node, const UnitStack& rate_units); + virtual void setParameters( + const AnyMap& node, const UnitStack& rate_units) override; //! Sets params to contain, in order, \f[ (A, B) \f] /** * @deprecated To be removed after Cantera 3.0; superseded by getFalloffCoeffs() */ - virtual void getParameters(double* params) const; + virtual void getParameters(double* params) const override; - virtual void getParameters(AnyMap& node) const; + virtual void getParameters(AnyMap& node) const override; protected: //! parameter a in the Tsang F_cent formulation. Dimensionless diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index ef9112dc2b..a890a67337 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -10,18 +10,15 @@ #include "cantera/base/AnyMap.h" #include "cantera/base/Units.h" -#include "Arrhenius.h" -#include "ChebyshevRate.h" -#include "Custom.h" -#include "Falloff.h" -#include "InterfaceRate.h" -#include "PlogRate.h" +#include "ReactionRate.h" namespace Cantera { class Kinetics; class ThirdBody; +class ArrheniusRate; // @todo Remove after Cantera 3.0 +class FalloffRate; // @todo Remove after Cantera 3.0 //! @defgroup reactionGroup Reactions and reaction rates @@ -33,29 +30,31 @@ class Reaction public: Reaction() {} Reaction(const Composition& reactants, const Composition& products, - shared_ptr rate={}); + shared_ptr rate, shared_ptr tbody=nullptr); + Reaction(const std::string& equation, + shared_ptr rate, shared_ptr tbody=nullptr); - //! Construct a Reaction and it's corresponding ReactionRate based on AnyMap (YAML) + //! Construct a Reaction and corresponding ReactionRate based on AnyMap (YAML) //! input. Reaction(const AnyMap& node, const Kinetics& kin); virtual ~Reaction() {} //! The reactant side of the chemical equation for this reaction - virtual std::string reactantString() const; + std::string reactantString() const; //! The product side of the chemical equation for this reaction - virtual std::string productString() const; + std::string productString() const; //! The chemical equation for this reaction std::string equation() const; //! Set the reactants and products based on the reaction equation. If a Kinetics //! object is provided, it is used to check that all reactants and products exist. - virtual void setEquation(const std::string& equation, const Kinetics* kin=0); + void setEquation(const std::string& equation, const Kinetics* kin=0); - //! The type of reaction - virtual std::string type() const; + //! The type of reaction, including reaction rate information + std::string type() const; //! Calculate the units of the rate constant. These are determined by the units //! of the standard concentration of the reactant species' phases and the phase @@ -74,11 +73,11 @@ class Reaction //! Ensure that the rate constant and other parameters for this reaction are valid. //! @since New in Cantera 3.0. - virtual void check(); + void check(); //! Ensure that the rate constant and other parameters for this reaction are valid. //! @deprecated To be removed after Cantera 3.0. Replaceable by check. - virtual void validate() { + void validate() { warn_deprecated("Reaction::validate", "Deprecated in Cantera 3.0 and to be removed thereafter; replaceable " "by Reaction::check."); @@ -87,7 +86,7 @@ class Reaction //! Perform validation checks that need access to a complete Kinetics objects, for // example to retrieve information about reactant / product species. - virtual void validate(Kinetics& kin) { + void validate(Kinetics& kin) { if (m_rate) { m_rate->validate(equation(), kin); } @@ -102,7 +101,7 @@ class Reaction AnyMap parameters(bool withInput=true) const; //! Set up reaction based on AnyMap *node* - virtual void setParameters(const AnyMap& node, const Kinetics& kin); + void setParameters(const AnyMap& node, const Kinetics& kin); //! Get validity flag of reaction bool valid() const { @@ -175,33 +174,37 @@ class Reaction //! Set reaction rate pointer void setRate(shared_ptr rate); - //! Get pointer to third-body + //! Get pointer to third-body handler shared_ptr thirdBody() { return m_third_body; } + //! Check whether reaction involves third body collider + //! @since New in Cantera 3.0. + bool usesThirdBody() const { + return bool(m_third_body); + } + protected: //! Store the parameters of a Reaction needed to reconstruct an identical //! object using the newReaction(AnyMap&, Kinetics&) function. Does not //! include user-defined fields available in the #input map. - virtual void getParameters(AnyMap& reactionNode) const; + void getParameters(AnyMap& reactionNode) const; //! Flag indicating whether reaction is set up correctly bool m_valid = true; - //! Helper function returning vector of undeclared third body species - //! and a boolean expression indicating whether the third body is specified. - //! @note The function is used by the checkSpecies method and only needed as long as - //! there is no unified approach to handle third body collision partners. - //! @param kin Kinetics object - virtual std::pair, bool> - undeclaredThirdBodies(const Kinetics& kin) const; + //! Flag indicating that serialization uses explicit type + bool m_explicit_rate = false; + + //! Flag indicating that object was instantiated from reactant/product compositions + bool m_from_composition = false; //! Reaction rate used by generic reactions shared_ptr m_rate; - //! Relative efficiencies of third-body species in enhancing the reaction - //! rate (if applicable) + //! Relative efficiencies of third-body species in enhancing the reaction rate + //! (if applicable) shared_ptr m_third_body; }; @@ -210,34 +213,71 @@ class Reaction class ThirdBody { public: - explicit ThirdBody(double default_efficiency=1.0); - + explicit ThirdBody() {}; + ThirdBody(const std::string& third_body); ThirdBody(const AnyMap& node); + //! @deprecated To be removed after Cantera 3.0; instantiate using string instead + ThirdBody(double default_efficiency); + + //! Name of the third body collider + //! @since New in Cantera 3.0 + std::string name() const { + return m_name; + } + + //! Set name of the third body collider + //! @since New in Cantera 3.0 + void setName(const std::string& third_body); + //! Set third-body efficiencies from AnyMap *node* + //! @deprecated To be removed after Cantera 3.0; renamed to setParameters void setEfficiencies(const AnyMap& node); + //! Set third-body efficiencies from AnyMap *node* + //! @since New in Cantera 3.0 + void setParameters(const AnyMap& node); + + //! Get third-body efficiencies from AnyMap *node* + //! @since New in Cantera 3.0 + void getParameters(AnyMap& node) const; + //! Get the third-body efficiency for species *k* double efficiency(const std::string& k) const; + //! Suffix representing the third body collider in reaction equation, for example + //! `+ M` or `(+M)` + //! @since New in Cantera 3.0 + std::string collider() const; + + //! Verify that all species involved in collision efficiencies are defined in the + //! Kinetics object. The function returns true if all species are found, and raises + //! an exception unless the Kinetics object is configured to skip undeclared + //! species, in which case false is returned. + //! @param rxn Reaction object + //! @param kin Kinetics object + //! @since New in Cantera 3.0 + bool checkSpecies(const Reaction& rxn, const Kinetics& kin) const; + //! Map of species to third body efficiency Composition efficiencies; - //! The default third body efficiency for species not listed in - //! #efficiencies. - double default_efficiency; - - //! Input explicitly specifies collision partner - bool specified_collision_partner; + //! The default third body efficiency for species not listed in #efficiencies. + double default_efficiency = 1.; //! Third body is used by law of mass action //! (`true` for three-body reactions, `false` for falloff reactions) - bool mass_action; + bool mass_action = true; + +protected: + //! Name of the third body collider + std::string m_name = "M"; }; //! A reaction with a non-reacting third body "M" that acts to add or remove //! energy from the reacting species +//! @deprecated To be removed after Cantera 3.0. Merged with Reaction class ThreeBodyReaction : public Reaction { public: @@ -246,18 +286,6 @@ class ThreeBodyReaction : public Reaction const ArrheniusRate& rate, const ThirdBody& tbody); ThreeBodyReaction(const AnyMap& node, const Kinetics& kin); - - virtual std::string type() const { - return "three-body"; - } - - virtual void setEquation(const std::string& equation, const Kinetics* kin=0); - bool detectEfficiencies(); - virtual void setParameters(const AnyMap& node, const Kinetics& kin); - virtual void getParameters(AnyMap& reactionNode) const; - - virtual std::string reactantString() const; - virtual std::string productString() const; }; @@ -267,23 +295,15 @@ class ThreeBodyReaction : public Reaction //! decreases as pressure increases due to collisional stabilization of a reaction //! intermediate; in this case, the forward rate constant is written as being //! proportional to the low-pressure rate constant. +//! @deprecated To be removed after Cantera 3.0. Merged with Reaction class FalloffReaction : public Reaction { public: FalloffReaction(); FalloffReaction(const Composition& reactants, const Composition& products, - const ReactionRate& rate, const ThirdBody& tbody); + const FalloffRate& rate, const ThirdBody& tbody); FalloffReaction(const AnyMap& node, const Kinetics& kin); - - virtual std::string type() const; - - virtual void setEquation(const std::string& equation, const Kinetics* kin); - virtual void setParameters(const AnyMap& node, const Kinetics& kin); - virtual void getParameters(AnyMap& reactionNode) const; - - virtual std::string reactantString() const; - virtual std::string productString() const; }; diff --git a/include/cantera/kinetics/ReactionRate.h b/include/cantera/kinetics/ReactionRate.h index 5e9e76623f..0619059221 100644 --- a/include/cantera/kinetics/ReactionRate.h +++ b/include/cantera/kinetics/ReactionRate.h @@ -85,6 +85,11 @@ class ReactionRate //! String identifying reaction rate specialization virtual const std::string type() const = 0; + //! String identifying sub-type of reaction rate specialization + virtual const std::string subType() const { + return ""; + } + //! Set parameters //! @param node AnyMap object containing reaction rate specification //! @param units unit definitions specific to rate information diff --git a/interfaces/cython/cantera/_utils.pxd b/interfaces/cython/cantera/_utils.pxd index a5e11da208..bb3abb8193 100644 --- a/interfaces/cython/cantera/_utils.pxd +++ b/interfaces/cython/cantera/_utils.pxd @@ -90,5 +90,5 @@ cdef pystr(string x) cdef comp_map_to_dict(Composition m) cdef Composition comp_map(X) except * -cdef CxxAnyMap dict_to_anymap(data) except * +cdef CxxAnyMap dict_to_anymap(dict data, cbool hyphenize=*) except * cdef anymap_to_dict(CxxAnyMap& m) diff --git a/interfaces/cython/cantera/_utils.pyx b/interfaces/cython/cantera/_utils.pyx index 5283c24c5b..4a402d0159 100644 --- a/interfaces/cython/cantera/_utils.pyx +++ b/interfaces/cython/cantera/_utils.pyx @@ -170,8 +170,18 @@ cdef anymap_to_dict(CxxAnyMap& m): return {pystr(item.first): anyvalue_to_python(item.first, item.second) for item in m.ordered()} -cdef CxxAnyMap dict_to_anymap(data) except *: +cdef CxxAnyMap dict_to_anymap(dict data, cbool hyphenize=False) except *: cdef CxxAnyMap m + if hyphenize: + # replace "_" by "-": while Python dictionaries typically use "_" in key names, + # the YAML convention uses "-" in field names + def _hyphenize(data): + if isinstance(data, dict): + return {k.replace("_", "-"): _hyphenize(v) for k, v in data.items()} + return data + + data = _hyphenize(data) + for k, v in data.items(): m[stringify(k)] = python_to_anyvalue(v, k) return m diff --git a/interfaces/cython/cantera/examples/kinetics/custom_reactions.py b/interfaces/cython/cantera/examples/kinetics/custom_reactions.py index 0970362786..51940fab0e 100644 --- a/interfaces/cython/cantera/examples/kinetics/custom_reactions.py +++ b/interfaces/cython/cantera/examples/kinetics/custom_reactions.py @@ -25,8 +25,7 @@ custom_reactions = [r for r in reactions] custom_reactions[2] = ct.Reaction( equation='H2 + O <=> H + OH', - rate=lambda T: 38.7 * T**2.7 * exp(-3150.15428/T), - kinetics=gas0) + rate=lambda T: 38.7 * T**2.7 * exp(-3150.15428/T)) gas1 = ct.Solution(thermo='ideal-gas', kinetics='gas', species=species, reactions=custom_reactions) diff --git a/interfaces/cython/cantera/reaction.pxd b/interfaces/cython/cantera/reaction.pxd index a956713480..0c0a836b50 100644 --- a/interfaces/cython/cantera/reaction.pxd +++ b/interfaces/cython/cantera/reaction.pxd @@ -31,6 +31,7 @@ cdef extern from "cantera/kinetics/ReactionRate.h" namespace "Cantera": cdef cppclass CxxReactionRate "Cantera::ReactionRate": CxxReactionRate() string type() + string subType() double eval(double) except +translate_exception double eval(double, double) except +translate_exception double eval(double, vector[double]&) except +translate_exception @@ -88,10 +89,49 @@ cdef extern from "cantera/cython/wrappers.h": cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxKinetics "Cantera::Kinetics" - cdef shared_ptr[CxxReaction] CxxNewReaction "Cantera::newReaction" (string) except +translate_exception cdef shared_ptr[CxxReaction] CxxNewReaction "newReaction" (CxxAnyMap&, CxxKinetics&) except +translate_exception cdef vector[shared_ptr[CxxReaction]] CxxGetReactions "getReactions" (CxxAnyValue&, CxxKinetics&) except +translate_exception + cdef cppclass CxxThirdBody "Cantera::ThirdBody": + CxxThirdBody() + CxxThirdBody(string) + double efficiency(string) + string name() + Composition efficiencies + double default_efficiency + cbool mass_action + + cdef cppclass CxxReaction "Cantera::Reaction": + CxxReaction() + CxxReaction(Composition&, Composition&, shared_ptr[CxxReactionRate]) except +translate_exception + CxxReaction(Composition&, Composition&, shared_ptr[CxxReactionRate], shared_ptr[CxxThirdBody]) except +translate_exception + CxxReaction(string&, shared_ptr[CxxReactionRate]) except +translate_exception + CxxReaction(string&, shared_ptr[CxxReactionRate], shared_ptr[CxxThirdBody]) except +translate_exception + + string reactantString() + string productString() + string equation() + void setEquation(const string&) except +translate_exception + string type() + CxxAnyMap parameters(cbool) except +translate_exception + cbool usesThirdBody() + CxxAnyMap input + Composition reactants + Composition products + Composition orders + string id + cbool reversible + cbool duplicate + cbool allow_nonreactant_orders + cbool allow_negative_orders + shared_ptr[CxxThirdBody] thirdBody() + CxxUnits rate_units + + shared_ptr[CxxReactionRate] rate() + void setRate(shared_ptr[CxxReactionRate]) + + +cdef extern from "cantera/kinetics/Falloff.h" namespace "Cantera": cdef cppclass CxxFalloffRate "Cantera::FalloffRate" (CxxReactionRate): CxxFalloffRate() CxxFalloffRate(CxxAnyMap) except +translate_exception @@ -119,12 +159,16 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxTsangRate "Cantera::TsangRate" (CxxFalloffRate): CxxTsangRate(CxxAnyMap) except +translate_exception + +cdef extern from "cantera/kinetics/PlogRate.h" namespace "Cantera": cdef cppclass CxxPlogRate "Cantera::PlogRate" (CxxReactionRate): CxxPlogRate() CxxPlogRate(CxxAnyMap) except +translate_exception CxxPlogRate(multimap[double, CxxArrheniusRate]) multimap[double, CxxArrheniusRate] getRates() + +cdef extern from "cantera/kinetics/ChebyshevRate.h" namespace "Cantera": cdef cppclass CxxChebyshevRate "Cantera::ChebyshevRate" (CxxReactionRate): CxxChebyshevRate() CxxChebyshevRate(CxxAnyMap) except +translate_exception @@ -137,10 +181,14 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": size_t nPressure() CxxArray2D& data() + +cdef extern from "cantera/kinetics/Custom.h" namespace "Cantera": cdef cppclass CxxCustomFunc1Rate "Cantera::CustomFunc1Rate" (CxxReactionRate): CxxCustomFunc1Rate() void setRateFunction(shared_ptr[CxxFunc1]) except +translate_exception + +cdef extern from "cantera/kinetics/InterfaceRate.h" namespace "Cantera": cdef cppclass CxxInterfaceRateBase "Cantera::InterfaceRateBase": void getCoverageDependencies(CxxAnyMap) void setCoverageDependencies(CxxAnyMap) except +translate_exception @@ -180,63 +228,6 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": CxxStickingBlowersMaselRate(CxxAnyMap) except +translate_exception CxxStickingBlowersMaselRate(double, double, double, double) - cdef cppclass CxxThirdBody "Cantera::ThirdBody": - CxxThirdBody() - CxxThirdBody(double) - double efficiency(string) - Composition efficiencies - double default_efficiency - - cdef cppclass CxxReaction "Cantera::Reaction": - CxxReaction() - - string reactantString() - string productString() - string equation() - void setEquation(const string&) except +translate_exception - string type() - CxxAnyMap parameters(cbool) except +translate_exception - CxxAnyMap input - Composition reactants - Composition products - Composition orders - string id - cbool reversible - cbool duplicate - cbool allow_nonreactant_orders - cbool allow_negative_orders - shared_ptr[CxxThirdBody] thirdBody() - CxxUnits rate_units - - shared_ptr[CxxReactionRate] rate() - void setRate(shared_ptr[CxxReactionRate]) - - cdef cppclass CxxFalloff "Cantera::FalloffRate": - CxxFalloff() - void updateTemp(double, double*) - double F(double, double*) - size_t workSize() - - size_t nParameters() - string type() - void getParameters(double*) - - cdef cppclass CxxChebyshev "Cantera::ChebyshevRate": - CxxChebyshev(double, double, double, double, CxxArray2D) - double Tmin() - double Tmax() - double Pmin() - double Pmax() - size_t nTemperature() - size_t nPressure() - CxxArray2D& data() - - cdef cppclass CxxThreeBodyReaction "Cantera::ThreeBodyReaction" (CxxReaction): - CxxThreeBodyReaction() - - cdef cppclass CxxFalloffReaction "Cantera::FalloffReaction" (CxxReaction): - CxxFalloffReaction() - cdef class ReactionRate: cdef shared_ptr[CxxReactionRate] _rate @@ -262,6 +253,12 @@ cdef class InterfaceRateBase(ArrheniusRateBase): cdef class StickRateBase(InterfaceRateBase): cdef CxxStickingCoverage* stick +cdef class ThirdBody: + cdef shared_ptr[CxxThirdBody] _third_body + cdef CxxThirdBody* third_body + @staticmethod + cdef wrap(shared_ptr[CxxThirdBody]) + cdef class Reaction: cdef shared_ptr[CxxReaction] _reaction cdef CxxReaction* reaction @@ -275,7 +272,3 @@ cdef class Arrhenius: cdef Reaction reaction # parent reaction, to prevent garbage collection @staticmethod cdef wrap(CxxArrheniusRate*) - -cdef class Falloff: - cdef shared_ptr[CxxFalloff] _falloff - cdef CxxFalloff* falloff diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index 7b77c9ea92..63344c7982 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -10,10 +10,6 @@ from .kinetics cimport Kinetics from ._utils cimport * from .units cimport * -# dictionary to store reaction classes -cdef dict _reaction_class_registry = {} - - # dictionary to store reaction rate classes cdef dict _reaction_rate_class_registry = {} @@ -41,6 +37,11 @@ cdef class ReactionRate: def __get__(self): return pystr(self.rate.type()) + property sub_type: + """ Get the C++ ReactionRate sub-type """ + def __get__(self): + return pystr(self.rate.subType()) + @staticmethod cdef wrap(shared_ptr[CxxReactionRate] rate): """ @@ -57,8 +58,10 @@ cdef class ReactionRate: # update global reaction class registry register_subclasses(ReactionRate) - # identify class - rate_type = pystr(rate.get().type()) + # identify class (either subType or type) + rate_type = pystr(rate.get().subType()) + if rate_type == "": + rate_type = pystr(rate.get().type()) cls = _reaction_rate_class_registry.get(rate_type, ReactionRate) # wrap C++ reaction rate @@ -72,10 +75,10 @@ cdef class ReactionRate: self.rate = self._rate.get() @classmethod - def from_dict(cls, data): + def from_dict(cls, data, hyphenize=True): """ Create a `ReactionRate` object from a dictionary corresponding to its YAML - representation. + representation. By default, underscores in keys are replaced by hyphens. An example for the creation of a `ReactionRate` from a dictionary is:: @@ -90,7 +93,7 @@ cdef class ReactionRate: f"Class method 'from_dict' was invoked from '{cls.__name__}' but " "should be called from base class 'ReactionRate'") - cdef CxxAnyMap any_map = dict_to_anymap(data) + cdef CxxAnyMap any_map = dict_to_anymap(data, hyphenize=hyphenize) cxx_rate = CxxNewReactionRate(any_map) return ReactionRate.wrap(cxx_rate) @@ -318,7 +321,9 @@ cdef class TwoTempPlasmaRate(ArrheniusRateBase): return self.rate.eval(temperature, elec_temp) def _from_dict(self, dict input_data): - self._rate.reset(new CxxTwoTempPlasmaRate(dict_to_anymap(input_data))) + self._rate.reset( + new CxxTwoTempPlasmaRate(dict_to_anymap(input_data, hyphenize=True)) + ) def _from_parameters(self, A, b, Ea_gas, Ea_electron): self._rate.reset(new CxxTwoTempPlasmaRate(A, b, Ea_gas, Ea_electron)) @@ -443,7 +448,9 @@ cdef class LindemannRate(FalloffRate): _reaction_rate_type = "Lindemann" def _from_dict(self, dict input_data): - self._rate.reset(new CxxLindemannRate(dict_to_anymap(input_data))) + self._rate.reset( + new CxxLindemannRate(dict_to_anymap(input_data, hyphenize=True)) + ) cdef set_cxx_object(self): self.rate = self._rate.get() @@ -461,7 +468,9 @@ cdef class TroeRate(FalloffRate): _reaction_rate_type = "Troe" def _from_dict(self, dict input_data): - self._rate.reset(new CxxTroeRate(dict_to_anymap(input_data))) + self._rate.reset( + new CxxTroeRate(dict_to_anymap(input_data, hyphenize=True)) + ) cdef set_cxx_object(self): self.rate = self._rate.get() @@ -479,7 +488,9 @@ cdef class SriRate(FalloffRate): _reaction_rate_type = "SRI" def _from_dict(self, dict input_data): - self._rate.reset(new CxxSriRate(dict_to_anymap(input_data))) + self._rate.reset( + new CxxSriRate(dict_to_anymap(input_data, hyphenize=True)) + ) cdef set_cxx_object(self): self.rate = self._rate.get() @@ -493,7 +504,9 @@ cdef class TsangRate(FalloffRate): _reaction_rate_type = "Tsang" def _from_dict(self, dict input_data): - self._rate.reset(new CxxTsangRate(dict_to_anymap(input_data))) + self._rate.reset( + new CxxTsangRate(dict_to_anymap(input_data, hyphenize=True)) + ) cdef set_cxx_object(self): self.rate = self._rate.get() @@ -992,6 +1005,87 @@ cdef class StickingBlowersMaselRate(StickRateBase): self.cxx_object().setDeltaH(delta_H) +cdef class ThirdBody: + r""" + Class representing third-body collision partners in three-body or falloff reactions. + + :param collider: + Name of the third-body collider. If ``M`` (default), the `default_efficiency` + is set to 1 and the collider is assumed to participate in a three-body reaction. + If the collider includes parentheses, - for example ``(+M)``, - a falloff form + is assumed, where the collider is not considered for the law of mass action. + For species other than ``M``, the third-body collider represents a named species + with collision efficiency 1, and the `default_efficiency` is set to zero. + :param efficiencies: + Non-default third-body efficiencies + :param default_efficiency: + Default third-body efficiency + + .. versionadded:: 3.0 + """ + def __cinit__(self, collider="M", *, + efficiencies=None, default_efficiency=None, init=True): + if not init: + return + self._third_body.reset(new CxxThirdBody(stringify(collider))) + self.third_body = self._third_body.get() + + if efficiencies is not None: + self.efficiencies = efficiencies + + if default_efficiency is not None: + self.default_efficiency = default_efficiency + + @staticmethod + cdef wrap(shared_ptr[CxxThirdBody] third_body): + tb = ThirdBody(init=False) + tb._third_body = third_body + tb.third_body = tb._third_body.get() + return tb + + property name: + """ + Get the name of the third-body collider used in the reaction equation. + """ + def __get__(self): + return pystr(self.third_body.name()) + + property mass_action: + """ + Retrieve flag indicating whether third-body collider participates + in the law of mass action. + """ + def __get__(self): + return self.third_body.mass_action + + property efficiencies: + """ + Get/Set a `dict` defining non-default third-body efficiencies for this reaction, + where the keys are the species names and the values are the efficiencies. + """ + def __get__(self): + return comp_map_to_dict(self.third_body.efficiencies) + def __set__(self, eff): + self.third_body.efficiencies = comp_map(eff) + + property default_efficiency: + """ + Get/Set the default third-body efficiency for this reaction, used for species + not in `efficiencies`. + """ + def __get__(self): + return self.third_body.default_efficiency + def __set__(self, default_eff): + self.third_body.default_efficiency = default_eff + + def efficiency(self, species): + """ + Get the efficiency of the third body named ``species`` considering both + the default efficiency and species-specific efficiencies. + """ + return self.third_body.efficiency(stringify(species)) + + cdef class Reaction: """ A class which stores data about a reaction and its rate parameterization so @@ -1049,73 +1143,108 @@ cdef class Reaction: rate-constant: {A: 3.87e+04 cm^3/mol/s, b: 2.7, Ea: 6260 cal/mol}}''', gas) """ - _reaction_type = "" def __cinit__(self, reactants=None, products=None, rate=None, *, - init=True, **kwargs): - if init: - rxn_type = self._reaction_type - self._reaction = CxxNewReaction(stringify((rxn_type))) - self.reaction = self._reaction.get() - if reactants: - self.reactants = reactants - if products: - self.products = products - - def __init__(self, reactants=None, products=None, rate=None, *, equation=None, - init=True, **kwargs): + equation=None, init=True, efficiencies=None, + Kinetics kinetics=None, third_body=None, **kwargs): + if kinetics: + warnings.warn( + "Parameter 'kinetics' is no longer used and will be removed after " + "Cantera 3.0.", DeprecationWarning) if not init: return - if equation: - self.reaction.setEquation(stringify(equation)) - - if isinstance(rate, dict): - if set(rate) == {"A", "b", "Ea"}: - # Allow simple syntax for Arrhenius rates - self.rate = ReactionRate.from_dict({"rate-constant": rate}) + cdef ReactionRate _rate + if isinstance(rate, ReactionRate): + _rate = rate + elif rate is None: + # default to Arrhenius expression + raise ValueError("Missing reaction rate information.") + elif isinstance(rate, dict): + if {"A", "b"} - set(rate) == set(): + # Allow simple syntax for Arrhenius-type rates + args = {"rate-constant": rate} + if "Ea0" in rate: + args.update({"type": "Blowers-Masel"}) + elif "Ea_gas" in rate: + args.update({"type": "two-temperature-plasma"}) + _rate = ReactionRate.from_dict(args) + else: + _rate = ReactionRate.from_dict(rate) + elif isinstance(rate, Arrhenius): + _rate = ArrheniusRate( + A=rate.pre_exponential_factor, + b=rate.temperature_exponent, + Ea=rate.activation_energy + ) + elif callable(rate): + _rate = CustomRate(rate) + else: + raise TypeError(f"Invalid rate definition with type '{type(rate)}'") + + cdef ThirdBody _third_body + if isinstance(third_body, ThirdBody): + _third_body = third_body + elif isinstance(third_body, str): + _third_body = ThirdBody(third_body) + elif efficiencies: + warnings.warn( + "Argument 'efficiencies' is deprecated and will be removed after " + "Cantera 3.0. Use ThirdBody instead.", DeprecationWarning) + third_body = "M" + _third_body = ThirdBody(third_body) + _third_body.efficiencies = efficiencies + + if reactants and products: + # create from reactant and product compositions + if third_body: + self._reaction.reset( + new CxxReaction(comp_map(reactants), comp_map(products), + _rate._rate, _third_body._third_body) + ) else: - self.rate = ReactionRate.from_dict(rate) - elif rate is not None: - self.rate = rate + self._reaction.reset( + new CxxReaction(comp_map(reactants), comp_map(products), + _rate._rate) + ) + elif equation: + # create from reaction equation + if third_body: + self._reaction.reset( + new CxxReaction(stringify(equation), + _rate._rate, _third_body._third_body) + ) + else: + self._reaction.reset( + new CxxReaction(stringify(equation), + _rate._rate) + ) + else: + # create default object + raise ValueError("Missing reactant and/or product information.") - if isinstance(self.rate, ReactionRate): - self.reaction.setRate(self._rate._rate) + self.reaction = self._reaction.get() + self._rate = _rate @staticmethod cdef wrap(shared_ptr[CxxReaction] reaction): """ Wrap a C++ Reaction object with a Python object of the correct derived type. """ - # ensure all reaction types are registered - if not _reaction_class_registry: - def register_subclasses(cls): - for c in cls.__subclasses__(): - rxn_type = getattr(c, "_reaction_type") - _reaction_class_registry[rxn_type] = c - register_subclasses(c) - - # update global reaction class registry - register_subclasses(Reaction) - - # identify class - rxn_type = pystr(reaction.get().type()) - cls = _reaction_class_registry.get(rxn_type, Reaction) - # wrap C++ reaction cdef Reaction R - R = cls(init=False) + R = Reaction(init=False) R._reaction = reaction R.reaction = R._reaction.get() R._rate = ReactionRate.wrap(R.reaction.rate()) return R @classmethod - def from_dict(cls, data, Kinetics kinetics): + def from_dict(cls, data, Kinetics kinetics, hyphenize=True): """ Create a `Reaction` object from a dictionary corresponding to its YAML - representation. + representation. By default, underscores in keys are replaced by hyphens. An example for the creation of a Reaction from a dictionary is:: @@ -1132,12 +1261,7 @@ cdef class Reaction: A `Kinetics` object whose associated phase(s) contain the species involved in the reaction. """ - if cls._reaction_type != "": - raise TypeError( - f"Class method 'from_dict' was invoked from '{cls.__name__}' but " - "should be called from base class 'Reaction'") - - cdef CxxAnyMap any_map = dict_to_anymap(data) + cdef CxxAnyMap any_map = dict_to_anymap(data, hyphenize=hyphenize) cxx_reaction = CxxNewReaction(any_map, deref(kinetics.kinetics)) return Reaction.wrap(cxx_reaction) @@ -1161,13 +1285,7 @@ cdef class Reaction: A `Kinetics` object whose associated phase(s) contain the species involved in the reaction. """ - if cls._reaction_type != "": - raise TypeError( - f"Class method 'from_yaml' was invoked from '{cls.__name__}' but " - "should be called from base class 'Reaction'") - - cdef CxxAnyMap any_map - any_map = AnyMapFromYamlString(stringify(text)) + cdef CxxAnyMap any_map = AnyMapFromYamlString(stringify(text)) cxx_reaction = CxxNewReaction(any_map, deref(kinetics.kinetics)) return Reaction.wrap(cxx_reaction) @@ -1357,7 +1475,7 @@ cdef class Reaction: self.reaction.input.clear() def __repr__(self): - return f"{self.equation} <{self.__class__.__name__}({self.rate.type})>" + return f"{self.equation} <{self.reaction_type}>" def __str__(self): return self.equation @@ -1368,6 +1486,84 @@ cdef class Reaction: cdef CxxUnits rate_units = self.reaction.rate_units return Units.copy(rate_units) + property third_body: + """ + Returns a `ThirdBody` object if `Reaction` uses a third body collider, and + ``None`` otherwise. + + .. versionadded:: 3.0 + """ + def __get__(self): + if self.reaction.usesThirdBody(): + return ThirdBody.wrap(self.reaction.thirdBody()) + + property efficiencies: + """ + Get/Set a `dict` defining non-default third-body efficiencies for this reaction, + where the keys are the species names and the values are the efficiencies. + + .. deprecated:: 3.0 + + To be removed after Cantera 3.0. Access via `third_body` property and + `ThirdBody` object instead. + """ + def __get__(self): + warnings.warn( + "Property 'efficiencies' is deprecated and will be removed after " + "Cantera 3.0. Use ThirdBody instead.", DeprecationWarning) + if self.third_body is None: + raise ValueError("Reaction does not involve third body collider") + return self.third_body.efficiencies + def __set__(self, eff): + warnings.warn( + "Property 'efficiencies' is deprecated and will be removed after " + "Cantera 3.0. Use ThirdBody instead.", DeprecationWarning) + if self.third_body is None: + raise ValueError("Reaction does not involve third body collider") + self.third_body.efficiencies = comp_map(eff) + + property default_efficiency: + """ + Get/Set the default third-body efficiency for this reaction, used for species + not in `efficiencies`. + + .. deprecated:: 3.0 + + To be removed after Cantera 3.0. Access via `third_body` property and + `ThirdBody` object instead. + """ + def __get__(self): + warnings.warn( + "Property 'default_efficiency' is deprecated and will be removed after " + "Cantera 3.0. Use ThirdBody instead.", DeprecationWarning) + if self.third_body is None: + raise ValueError("Reaction does not involve third body collider") + return self.third_body.default_efficiency + def __set__(self, default_eff): + warnings.warn( + "Property 'default_efficiency' is deprecated and will be removed after " + "Cantera 3.0. Use ThirdBody instead.", DeprecationWarning) + if self.third_body is None: + raise ValueError("Reaction does not involve third body collider") + self.third_body.default_efficiency = default_eff + + def efficiency(self, species): + """ + Get the efficiency of the third body named ``species`` considering both + the default efficiency and species-specific efficiencies. + + .. deprecated:: 3.0 + + To be removed after Cantera 3.0. Access via `third_body` property and + `ThirdBody` object instead. + """ + warnings.warn( + "Method 'efficiency' is deprecated and will be removed after " + "Cantera 3.0. Use ThirdBody instead.", DeprecationWarning) + if self.third_body is None: + raise ValueError("Reaction does not involve third body collider") + return self.third_body.efficiency(stringify(species)) + cdef class Arrhenius: r""" @@ -1462,86 +1658,15 @@ cdef class ThreeBodyReaction(Reaction): type: three-body rate-constant: {A: 1.2e+17 cm^6/mol^2/s, b: -1.0, Ea: 0.0 cal/mol} efficiencies: {H2: 2.4, H2O: 15.4, AR: 0.83} - """ - _reaction_type = "three-body" - - cdef CxxThreeBodyReaction* cxx_threebody(self): - return self.reaction - - cdef CxxThirdBody* thirdbody(self): - return (self.cxx_threebody().thirdBody().get()) - - def __init__(self, reactants=None, products=None, rate=None, *, equation=None, - efficiencies=None, Kinetics kinetics=None, init=True, **kwargs): - - if reactants and products and not equation: - equation = self.equation - - if isinstance(rate, ArrheniusRate): - self._reaction.reset(new CxxThreeBodyReaction()) - self.reaction = self._reaction.get() - if reactants and products: - self.reactants = reactants - self.products = products - else: - self.reaction.setEquation(stringify(equation)) - if efficiencies: - self.efficiencies = efficiencies - self.rate = rate - return - - if init and equation and kinetics: - rxn_type = self._reaction_type - spec = {"equation": equation, "type": rxn_type} - if isinstance(rate, dict): - spec["rate-constant"] = rate - elif isinstance(rate, Arrhenius) or rate is None: - spec["rate-constant"] = dict.fromkeys(["A", "b", "Ea"], 0.) - elif rate is None: - pass - elif isinstance(rate, ArrheniusRate): - pass - else: - raise TypeError("Invalid rate definition") - - if isinstance(efficiencies, dict): - spec["efficiencies"] = efficiencies - self._reaction = CxxNewReaction(dict_to_anymap(spec), - deref(kinetics.kinetics)) - self.reaction = self._reaction.get() - - self._rate = ReactionRate.wrap(self.reaction.rate()) - if isinstance(rate, Arrhenius): - self.rate = rate - - property efficiencies: - """ - Get/Set a `dict` defining non-default third-body efficiencies for this - reaction, where the keys are the species names and the values are the - efficiencies. - """ - def __get__(self): - return comp_map_to_dict(self.thirdbody().efficiencies) - def __set__(self, eff): - self.thirdbody().efficiencies = comp_map(eff) + .. deprecated:: 3.0 - property default_efficiency: - """ - Get/Set the default third-body efficiency for this reaction, used for - species used for species not in `efficiencies`. - """ - def __get__(self): - return self.thirdbody().default_efficiency - def __set__(self, default_eff): - self.thirdbody().default_efficiency = default_eff + Class to be removed after Cantera 3.0. Absorbed by `Reaction`. + """ - def efficiency(self, species): - """ - Get the efficiency of the third body named ``species`` considering both - the default efficiency and species-specific efficiencies. - """ - return self.thirdbody().efficiency(stringify(species)) + def __init__(self, *args, **kwargs): + warnings.warn("Class to be removed after Cantera 3.0; no specialization " + "necessary.", DeprecationWarning) cdef class FalloffReaction(Reaction): @@ -1567,83 +1692,15 @@ cdef class FalloffReaction(Reaction): high-P-rate-constant: {A: 7.4e+10, b: -0.37, Ea: 0.0 cal/mol} Troe: {A: 0.7346, T3: 94.0, T1: 1756.0, T2: 5182.0} efficiencies: {AR: 0.7, H2: 2.0, H2O: 6.0} - """ - _reaction_type = "falloff" - - cdef CxxFalloffReaction* cxx_object(self): - return self.reaction - - cdef CxxThirdBody* thirdbody(self): - return (self.cxx_object().thirdBody().get()) - - def __init__(self, reactants=None, products=None, rate=None, *, equation=None, - efficiencies=None, Kinetics kinetics=None, init=True, **kwargs): - - if reactants and products and not equation: - equation = self.equation - - if isinstance(rate, FalloffRate): - self._reaction.reset(new CxxFalloffReaction()) - self.reaction = self._reaction.get() - if reactants and products: - self.reactants = reactants - self.products = products - else: - self.reaction.setEquation(stringify(equation)) - if efficiencies: - self.efficiencies = efficiencies - self.rate = rate - return - - if init and equation and kinetics: - rxn_type = self._reaction_type - spec = {"equation": equation, "type": rxn_type} - if isinstance(rate, dict): - for key, value in rate.items(): - spec[key.replace("_", "-")] = value - elif isinstance(rate, FalloffRate) or rate is None: - pass - else: - raise TypeError( - f"Invalid rate definition; type is '{type(rate).__name__}'") - - if isinstance(efficiencies, dict): - spec["efficiencies"] = efficiencies - - self._reaction = CxxNewReaction(dict_to_anymap(spec), - deref(kinetics.kinetics)) - self.reaction = self._reaction.get() - self._rate = ReactionRate.wrap(self.reaction.rate()) - - property efficiencies: - """ - Get/Set a `dict` defining non-default third-body efficiencies for this - reaction, where the keys are the species names and the values are the - efficiencies. - """ - def __get__(self): - return comp_map_to_dict(self.thirdbody().efficiencies) - def __set__(self, eff): - self.thirdbody().efficiencies = comp_map(eff) - - property default_efficiency: - """ - Get/Set the default third-body efficiency for this reaction, used for - species used for species not in `efficiencies`. - """ - def __get__(self): - return self.thirdbody().default_efficiency - def __set__(self, default_eff): - self.thirdbody().default_efficiency = default_eff + .. deprecated:: 3.0 - def efficiency(self, species): - """ - Get the efficiency of the third body named ``species`` considering both - the default efficiency and species-specific efficiencies. - """ - return self.thirdbody().efficiency(stringify(species)) + Class to be removed after Cantera 3.0. Absorbed by `Reaction`. + """ + def __init__(self, *args, **kwargs): + warnings.warn("Class to be removed after Cantera 3.0; no specialization " + "necessary.", DeprecationWarning) cdef class ChemicallyActivatedReaction(FalloffReaction): """ @@ -1652,8 +1709,10 @@ cdef class ChemicallyActivatedReaction(FalloffReaction): that the forward rate constant is written as being proportional to the low- pressure rate constant. """ - _reaction_type = "chemically-activated" + def __init__(self, *args, **kwargs): + warnings.warn("Class to be removed after Cantera 3.0; no specialization " + "necessary.", DeprecationWarning) cdef class CustomReaction(Reaction): """ @@ -1674,4 +1733,3 @@ cdef class CustomReaction(Reaction): def __init__(self, *args, **kwargs): warnings.warn("Class to be removed after Cantera 3.0; no specialization " "necessary.", DeprecationWarning) - super().__init__(*args, **kwargs) diff --git a/interfaces/cython/cantera/test/test_jacobian.py b/interfaces/cython/cantera/test/test_jacobian.py index 2f858a9714..8d17d546d4 100644 --- a/interfaces/cython/cantera/test/test_jacobian.py +++ b/interfaces/cython/cantera/test/test_jacobian.py @@ -386,7 +386,7 @@ class TestFalloff(HydrogenOxygen, utilities.CanteraTest): # Fall-off reaction rxn_idx = 21 equation = "2 OH (+M) <=> H2O2 (+M)" - rate_type = "Troe" + rate_type = "falloff" rtol = 1e-4 diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index c428fe988a..b0744157de 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -63,17 +63,17 @@ def test_legacy_reaction_rate(self): fwd_rates_legacy = self.phase.forward_rate_constants ct.use_legacy_rate_constants(False) fwd_rates = self.phase.forward_rate_constants - ix_3b = np.array([r.reaction_type == "three-body" for r in self.phase.reactions()]) + ix_3b = np.array([r.reaction_type == "three-body-Arrhenius" for r in self.phase.reactions()]) ix_other = ix_3b == False self.assertArrayNear(fwd_rates_legacy[ix_other], fwd_rates[ix_other]) self.assertFalse((fwd_rates_legacy[ix_3b] == fwd_rates[ix_3b]).any()) def test_reaction_type(self): - self.assertIn(self.phase.reaction(0).reaction_type, "three-body") - self.assertIn(self.phase.reaction(2).reaction_type, "reaction") + self.assertIn(self.phase.reaction(0).reaction_type, "three-body-Arrhenius") + self.assertIn(self.phase.reaction(2).reaction_type, "Arrhenius") self.assertIn(self.phase.reaction(2).rate.type, "Arrhenius") - self.assertEqual(self.phase.reaction(21).reaction_type, "falloff") + self.assertEqual(self.phase.reaction(21).reaction_type, "falloff-Troe") with self.assertRaisesRegex(ct.CanteraError, "outside valid range"): self.phase.reaction(33).reaction_type @@ -1115,10 +1115,10 @@ def test_from_yaml(self): " efficiencies: {H2: 2.4, H2O: 15.4, AR: 0.83}}", self.gas) - self.assertTrue(isinstance(r, ct.ThreeBodyReaction)) + assert r.third_body is not None self.assertEqual(r.reactants['O'], 2) self.assertEqual(r.products['O2'], 1) - self.assertEqual(r.efficiencies['H2O'], 15.4) + self.assertEqual(r.third_body.efficiencies['H2O'], 15.4) self.assertEqual(r.rate.temperature_exponent, -1.0) self.assertIn('O', r) self.assertIn('O2', r) @@ -1197,7 +1197,6 @@ def test_negative_A(self): def test_negative_A_falloff(self): species = ct.Species.list_from_file("gri30.yaml") - r = ct.FalloffReaction('NH:1, NO:1', 'N2O:1, H:1') low_rate = ct.Arrhenius(2.16e13, -0.23, 0) high_rate = ct.Arrhenius(-8.16e12, -0.5, 0) @@ -1213,15 +1212,15 @@ def test_negative_A_falloff(self): rate.low_rate = low_rate rate.low_rate = ct.Arrhenius(-2.16e13, -0.23, 0) - r.rate = rate + rxn = ct.Reaction("NH:1, NO:1", "N2O:1, H:1", rate) gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', - species=species, reactions=[r]) + species=species, reactions=[rxn]) self.assertLess(gas.forward_rate_constants, 0) def test_threebody(self): - r = ct.ThreeBodyReaction({"O":1, "H":1}, {"OH":1}, - ct.ArrheniusRate(5e11, -1.0, 0.0)) - r.efficiencies = {"AR":0.7, "H2":2.0, "H2O":6.0} + tb = ct.ThirdBody(efficiencies={"AR":0.7, "H2":2.0, "H2O":6.0}) + r = ct.Reaction({"O":1, "H":1}, {"OH":1}, + ct.ArrheniusRate(5e11, -1.0, 0.0), third_body=tb) gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=self.species, reactions=[r]) @@ -1235,10 +1234,11 @@ def test_threebody(self): def test_falloff(self): high_rate = ct.Arrhenius(7.4e10, -0.37, 0.0) low_rate = ct.Arrhenius(2.3e12, -0.9, -1700 * 1000 * 4.184) - r = ct.FalloffReaction("OH:2", "H2O2:1", - ct.TroeRate(low_rate, high_rate, [0.7346, 94, 1756, 5182])) - r.efficiencies = {"AR":0.7, "H2":2.0, "H2O":6.0} - self.assertEqual(r.rate.type, "Troe") + tb = ct.ThirdBody(efficiencies={"AR":0.7, "H2":2.0, "H2O":6.0}) + r = ct.Reaction("OH:2", "H2O2:1", + ct.TroeRate(low_rate, high_rate, [0.7346, 94, 1756, 5182]), + third_body=tb) + self.assertEqual(r.rate.type, "falloff") gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=self.species, reactions=[r]) @@ -1561,12 +1561,12 @@ def test_modify_falloff(self): gas = ct.Solution('gri30.yaml', transport_model=None) gas.TPX = 1100, 3 * ct.one_atm, 'CH4:1.0, O2:0.4, CO2:0.1, H2O:0.05' r0 = gas.reaction(11) - self.assertEqual(r0.rate.type, "Lindemann") + self.assertEqual(r0.rate.type, "falloff") # these two reactions happen to have the same third-body efficiencies r1 = gas.reaction(49) r2 = gas.reaction(53) - self.assertEqual(r2.rate.type, "Troe") - self.assertEqual(r1.efficiencies, r2.efficiencies) + self.assertEqual(r2.rate.type, "falloff") + self.assertEqual(r1.third_body.efficiencies, r2.third_body.efficiencies) r2.rate = r1.rate gas.modify_reaction(53, r2) diff --git a/interfaces/cython/cantera/test/test_reaction.py b/interfaces/cython/cantera/test/test_reaction.py index 3ffdf60e2b..c05896c6fc 100644 --- a/interfaces/cython/cantera/test/test_reaction.py +++ b/interfaces/cython/cantera/test/test_reaction.py @@ -24,20 +24,19 @@ def test_implicit_three_body(self): rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} """ rxn1 = ct.Reaction.from_yaml(yaml1, self.gas) - self.assertEqual(rxn1.reaction_type, "three-body") - self.assertEqual(rxn1.default_efficiency, 0.) - self.assertEqual(rxn1.efficiencies, {"O2": 1}) + self.assertEqual(rxn1.reaction_type, "three-body-Arrhenius") + self.assertEqual(rxn1.third_body.default_efficiency, 0.) + self.assertEqual(rxn1.third_body.efficiencies, {"O2": 1}) yaml2 = """ equation: H + O2 + M <=> HO2 + M rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} - type: three-body default-efficiency: 0 efficiencies: {O2: 1.0} """ rxn2 = ct.Reaction.from_yaml(yaml2, self.gas) - self.assertEqual(rxn1.efficiencies, rxn2.efficiencies) - self.assertEqual(rxn1.default_efficiency, rxn2.default_efficiency) + self.assertEqual(rxn1.third_body.efficiencies, rxn2.third_body.efficiencies) + self.assertEqual(rxn1.third_body.default_efficiency, rxn2.third_body.default_efficiency) def test_duplicate(self): # @todo simplify this test @@ -54,7 +53,6 @@ def test_duplicate(self): yaml2 = """ equation: H + O2 + M <=> HO2 + M rate-constant: {A: 1.126e+19, b: -0.76, Ea: 0.0} - type: three-body default-efficiency: 0 efficiencies: {H2O: 1} """ @@ -63,8 +61,8 @@ def test_duplicate(self): self.assertEqual(rxn1.reaction_type, rxn2.reaction_type) self.assertEqual(rxn1.reactants, rxn2.reactants) self.assertEqual(rxn1.products, rxn2.products) - self.assertEqual(rxn1.efficiencies, rxn2.efficiencies) - self.assertEqual(rxn1.default_efficiency, rxn2.default_efficiency) + self.assertEqual(rxn1.third_body.efficiencies, rxn2.third_body.efficiencies) + self.assertEqual(rxn1.third_body.default_efficiency, rxn2.third_body.default_efficiency) gas1.add_reaction(rxn1) gas1.add_reaction(rxn2) @@ -97,7 +95,7 @@ def test_non_integer_stoich(self): rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} """ rxn = ct.Reaction.from_yaml(yaml, self.gas) - self.assertEqual(rxn.reaction_type, "reaction") + self.assertEqual(rxn.reaction_type, "Arrhenius") def test_not_three_body(self): # check that insufficient reactants prevent automatic conversion @@ -106,7 +104,7 @@ def test_not_three_body(self): rate-constant: {A: 2.1e+15, b: -0.69, Ea: 2850.0} """ rxn = ct.Reaction.from_yaml(yaml, self.gas) - self.assertEqual(rxn.reaction_type, "reaction") + self.assertEqual(rxn.reaction_type, "Arrhenius") def test_user_override(self): # check that type specification prevents automatic conversion @@ -116,7 +114,9 @@ def test_user_override(self): type: elementary """ rxn = ct.Reaction.from_yaml(yaml, self.gas) - self.assertEqual(rxn.reaction_type, "reaction") + self.assertEqual(rxn.reaction_type, "Arrhenius") + assert "type" in rxn.input_data + assert rxn.input_data["type"] == "elementary" class ReactionRateTests: @@ -409,6 +409,7 @@ def test_from_parts(self): class FalloffRateTests(ReactionRateTests): # test Falloff rate expressions + _type = "falloff" _n_data = [0] # list of valid falloff coefficient array lengths @classmethod @@ -469,7 +470,6 @@ class TestLindemannRate(FalloffRateTests, utilities.CanteraTest): # test Lindemann rate expressions _cls = ct.LindemannRate - _type = "Lindemann" _index = 7 _parts = { "falloff_coeffs": [], @@ -495,7 +495,6 @@ class TestTroeRate(FalloffRateTests, utilities.CanteraTest): # test Troe rate expressions _cls = ct.TroeRate - _type = "Troe" _index = 2 _parts = {"falloff_coeffs": [0.7346, 94.0, 1756.0, 5182.0]} _input = { @@ -528,7 +527,6 @@ class TestSriRate(FalloffRateTests, utilities.CanteraTest): # test SRI rate expressions _cls = ct.SriRate - _type = "SRI" _index = 8 _parts = {"falloff_coeffs": [1.1, 700.0, 1234.0, 56.0, 0.7]} _input = { @@ -550,7 +548,6 @@ class TestTsangRate(FalloffRateTests, utilities.CanteraTest): # test Tsang rate expressions _cls = ct.TsangRate - _type = "Tsang" _index = 9 _parts = {"falloff_coeffs": [0.95, -1.0e-04]} _input = { @@ -935,10 +932,10 @@ class ReactionTests: _rate_cls = None # corresponding reaction rate type _equation = None # reaction equation string _rate = None # parameters for reaction rate object constructor + _3rd_body = None # object representing third-body collider _rate_obj = None # reaction rate object _kwargs = {} # additional parameters required by constructor _index = None # index of reaction in "kineticsfromscratch.yaml" - _rxn_type = "reaction" # name of reaction type _rate_type = None # name of reaction rate type _yaml = None # YAML parameterization @@ -973,29 +970,22 @@ def from_dict(self): return self.finalize(rxn) 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() - rxn = self._cls(equation=self._equation, rate=rate, kinetics=self.soln, - **self._kwargs) + # create reaction object with an "empty" rate of the correct type + rxn = ct.Reaction(equation=self._equation, + rate=self._rate_cls(), third_body=self._3rd_body) return self.finalize(rxn) def from_rate(self, rate): - 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") - rxn = self._cls(equation=self._equation, rate=rate, kinetics=self.soln, - **self._kwargs) + # create reaction object from dictionary + rxn = ct.Reaction(equation=self._equation, + rate=rate, third_body=self._3rd_body) return self.finalize(rxn) def from_parts(self): # create reaction rate object from parts orig = self.soln.reaction(self._index) - rxn = self._cls(orig.reactants, orig.products, rate=self._rate_obj) + rxn = ct.Reaction(orig.reactants, orig.products, + rate=self._rate_obj, third_body=self._3rd_body) rxn.reversible = "<=>" in self._equation return self.finalize(rxn) @@ -1086,8 +1076,6 @@ def test_raises_invalid_rate(self): def test_no_rate(self): # check behavior for instantiation from keywords / no rate - if self._rate_obj is None: - return rxn = self.from_empty() self.assertIsNaN(self.eval_rate(rxn.rate)) @@ -1128,12 +1116,10 @@ def check_equal(self, one, two): class TestElementary(ReactionTests, utilities.CanteraTest): # test elementary reaction - _cls = ct.Reaction _rate_cls = ct.ArrheniusRate _equation = "H2 + O <=> H + OH" _rate = {"A": 38.7, "b": 2.7, "Ea": 2.619184e+07} _index = 0 - _rxn_type = "reaction" _rate_type = "Arrhenius" _yaml = """ equation: O + H2 <=> H + OH @@ -1150,30 +1136,29 @@ def setUpClass(cls): class TestThreeBody(TestElementary): # test three-body reaction - _cls = ct.ThreeBodyReaction _equation = "2 O + M <=> O2 + M" _rate = {"A": 1.2e11, "b": -1.0, "Ea": 0.0} - _kwargs = {"efficiencies": {"H2": 2.4, "H2O": 15.4, "AR": 0.83}} + _3rd_body = ct.ThirdBody(efficiencies={"H2": 2.4, "H2O": 15.4, "AR": 0.83}) _index = 1 - _rxn_type = "three-body" _yaml = """ equation: 2 O + M <=> O2 + M - type: three-body + # type: three-body # optional rate-constant: {A: 1.2e+11, b: -1.0, Ea: 0.0 cal/mol} efficiencies: {H2: 2.4, H2O: 15.4, AR: 0.83} """ - def from_parts(self): - rxn = ReactionTests.from_parts(self) - rxn.efficiencies = self._kwargs["efficiencies"] - return rxn - def test_efficiencies(self): # check efficiencies - rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.soln, - **self._kwargs) + rxn = self.from_rate(self._rate_obj) + self.assertEqual(rxn.third_body.efficiencies, self._3rd_body.efficiencies) - self.assertEqual(rxn.efficiencies, self._kwargs["efficiencies"]) + def test_serialization_type(self): + # test serialization output + rxn = self.from_yaml() + assert "type" not in rxn.input_data + orig = self.soln.reaction(self._index) + assert "type" in orig.input_data + assert orig.input_data["type"] == "three-body" class TestImplicitThreeBody(TestThreeBody): @@ -1181,24 +1166,27 @@ class TestImplicitThreeBody(TestThreeBody): _equation = "H + 2 O2 <=> HO2 + O2" _rate = {"A": 2.08e+19, "b": -1.24, "Ea": 0.0} - _kwargs = {} + _3rd_body = ct.ThirdBody("O2") _index = 5 _yaml = """ equation: H + 2 O2 <=> HO2 + O2 rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} """ - def from_parts(self): - rxn = ReactionTests.from_parts(self) - rxn.efficiencies = {"O2": 1.} - rxn.default_efficiency = 0 - return rxn - def test_efficiencies(self): # overload of default tester rxn = self.from_rate(self._rate_obj) - self.assertEqual(rxn.efficiencies, {"O2": 1.}) - self.assertEqual(rxn.default_efficiency, 0.) + self.assertEqual(rxn.third_body.efficiencies, {"O2": 1.}) + self.assertEqual(rxn.third_body.default_efficiency, 0.) + + def test_serialization_type(self): + # test serialization output + rxn = self.from_yaml() + assert "type" not in rxn.input_data + assert "efficiencies" not in rxn.input_data + orig = self.soln.reaction(self._index) + assert "type" not in orig.input_data + assert "efficiencies" not in orig.input_data class TestTwoTempPlasma(ReactionTests, utilities.CanteraTest): @@ -1207,6 +1195,7 @@ class TestTwoTempPlasma(ReactionTests, utilities.CanteraTest): _rate_cls = ct.TwoTempPlasmaRate _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 = """ @@ -1226,11 +1215,12 @@ def test_reversible(self): class TestBlowersMasel(ReactionTests, utilities.CanteraTest): - # test updated version of Blowers-Masel reaction + # test elementary version of Blowers-Masel reaction _rate_cls = ct.BlowersMaselRate _rate_type = "Blowers-Masel" _equation = "O + H2 <=> H + OH" + _rate = {"A": 38700, "b": 2.7, "Ea0": 1.0958665856e8, "w": 1.7505856e13} _rate_obj = ct.BlowersMaselRate(A=38700, b=2.7, Ea0=1.0958665856e8, w=1.7505856e13) _index = 6 _yaml = """ @@ -1245,10 +1235,23 @@ def eval_rate(self, rate): return rate(self.soln.T) +class TestThreeBodyBlowersMasel(TestBlowersMasel): + # test three-body version of Blowers-Masel reaction + + _equation = "O + H2 + M <=> H2O + M" + _3rd_body = ct.ThirdBody() + _index = 13 + _yaml = """ + equation: O + H2 + M <=> H2O + M + type: Blowers-Masel + rate-constant: {A: 38700, b: 2.7, Ea0: 2.619184e4 cal/mol, w: 4.184e9 cal/mol} + """ + + class TestTroe(ReactionTests, utilities.CanteraTest): # test Troe falloff reaction - _cls = ct.FalloffReaction + _rate_cls = ct.TroeRate _equation = "2 OH (+ M) <=> H2O2 (+ M)" _rate = { "type": "falloff", @@ -1256,10 +1259,9 @@ class TestTroe(ReactionTests, utilities.CanteraTest): "high_P_rate_constant": {"A": 7.4e+10, "b": -0.37, "Ea": 0.0}, "Troe": {"A": 0.7346, "T3": 94.0, "T1": 1756.0, "T2": 5182.0} } - _kwargs = {"efficiencies": {"AR": 0.7, "H2": 2.0, "H2O": 6.0}} + _3rd_body = ct.ThirdBody("(+M)", efficiencies={"AR": 0.7, "H2": 2.0, "H2O": 6.0}) _index = 2 - _rxn_type = "falloff" - _rate_type = "Troe" + _rate_type = "falloff" _yaml = """ equation: 2 OH (+ M) <=> H2O2 (+ M) # Reaction 3 type: falloff @@ -1284,26 +1286,20 @@ def eval_rate(self, rate): concm = self.soln.third_body_concentrations[self._index] return rate(self.soln.T, concm) - def from_parts(self): - rxn = ReactionTests.from_parts(self) - rxn.efficiencies = self._kwargs["efficiencies"] - return rxn - class TestLindemann(ReactionTests, utilities.CanteraTest): # test Lindemann falloff reaction - _cls = ct.FalloffReaction + _rate_cls = ct.LindemannRate _equation = "2 OH (+ M) <=> H2O2 (+ M)" _rate = { "type": "falloff", "low_P_rate_constant": {"A": 2.3e+12, "b": -0.9, "Ea": -7112800.0}, "high_P_rate_constant": {"A": 7.4e+10, "b": -0.37, "Ea": 0.0} } - _kwargs = {"efficiencies": {"AR": 0.7, "H2": 2.0, "H2O": 6.0}} + _3rd_body = ct.ThirdBody("(+M)", efficiencies={"AR": 0.7, "H2": 2.0, "H2O": 6.0}) _index = 7 - _rxn_type = "falloff" - _rate_type = "Lindemann" + _rate_type = "falloff" _yaml = """ equation: 2 OH (+ M) <=> H2O2 (+ M) # Reaction 8 duplicate: true @@ -1326,16 +1322,11 @@ def eval_rate(self, rate): concm = self.soln.third_body_concentrations[self._index] return rate(self.soln.T, concm) - def from_parts(self): - rxn = ReactionTests.from_parts(self) - rxn.efficiencies = self._kwargs["efficiencies"] - return rxn - class TestChemicallyActivated(ReactionTests, utilities.CanteraTest): # test Chemically Activated falloff reaction - _cls = ct.FalloffReaction + _rate_cls = ct.LindemannRate _equation = "H2O + OH (+M) <=> HO2 + H2 (+M)" _rate = { "type": "chemically-activated", @@ -1343,8 +1334,7 @@ class TestChemicallyActivated(ReactionTests, utilities.CanteraTest): "high_P_rate_constant": {"A": 5.88E-14, "b": 6.721, "Ea": -12644997.768} } _index = 10 - _rxn_type = "chemically-activated" - _rate_type = "Lindemann" + _rate_type = "chemically-activated" _yaml = """ equation: H2O + OH (+M) <=> HO2 + H2 (+M) # Reaction 11 units: {length: cm, quantity: mol, activation-energy: cal/mol} @@ -1371,9 +1361,7 @@ def eval_rate(self, rate): class TestPlog(ReactionTests, utilities.CanteraTest): # test Plog reaction - _cls = ct.Reaction _rate_cls = ct.PlogRate - _rxn_type = "reaction" _rate_type = "pressure-dependent-Arrhenius" _equation = "H2 + O2 <=> 2 OH" _rate = { @@ -1416,12 +1404,17 @@ def check_rates(self, rates, other): class TestChebyshev(ReactionTests, utilities.CanteraTest): # test Chebyshev reaction - _cls = ct.Reaction _rate_cls = ct.ChebyshevRate - _rxn_type = "reaction" _rate_type = "Chebyshev" _equation = "HO2 <=> OH + O" - _rate = None + _rate = { + "type": "Chebyshev", + "temperature-range": [290., 3000.], + "pressure-range": [1000., 10000000.0], + "data": + [[8.2883, -1.1397, -0.12059, 0.016034], + [1.9764, 1.0037, 7.2865e-03, -0.030432], + [0.3177, 0.26889, 0.094806, -7.6385e-03]]} _rate_obj = ct.ChebyshevRate( temperature_range=(290., 3000.), pressure_range=(1000., 10000000.0), data=[[ 8.2883e+00, -1.1397e+00, -1.2059e-01, 1.6034e-02], @@ -1447,12 +1440,10 @@ class TestCustom(ReactionTests, utilities.CanteraTest): # test Custom reaction # probe O + H2 <=> H + OH - _cls = ct.Reaction _rate_cls = ct.CustomRate _equation = "H2 + O <=> H + OH" _rate_obj = ct.CustomRate(lambda T: 38.7 * T**2.7 * exp(-3150.15428/T)) _index = 0 - _rxn_type = "reaction" _rate_type = "custom-rate-function" _yaml = None diff --git a/interfaces/cython/cantera/yaml2ck.py b/interfaces/cython/cantera/yaml2ck.py index e4f4825e3c..81a8797a6d 100644 --- a/interfaces/cython/cantera/yaml2ck.py +++ b/interfaces/cython/cantera/yaml2ck.py @@ -367,6 +367,7 @@ def build_reactions_text(reactions: Iterable[ct.Reaction]): ) reaction_order += sum(reac.orders.values()) unit_conversion_factor = 1_000.0 ** (reaction_order - 1) + if reac.rate.type == "Chebyshev": reaction_lines.append( arrhenius_line.format( @@ -391,13 +392,11 @@ def build_reactions_text(reactions: Iterable[ct.Reaction]): coeffs[0, 0] += math.log10(unit_conversion_factor) for row in coeffs: reaction_lines.append(f"CHEB /{' '.join(map(str, row))}/") - elif reac.reaction_type in ("reaction", "three-body"): - if reac.reaction_type == "three-body": + + elif reac.rate.type == "Arrhenius": + if reac.reaction_type.startswith("three-body"): unit_conversion_factor *= 1_000.0 - if reac.rate.type == "pressure-dependent-Arrhenius": - rate = reac.rate.rates[-1][1] - else: - rate = reac.rate + rate = reac.rate reaction_lines.append( arrhenius_line.format( equation=reac.equation, @@ -407,85 +406,83 @@ def build_reactions_text(reactions: Iterable[ct.Reaction]): E_a=rate.activation_energy / CALORIES_CONSTANT, ) ) - if reac.rate.type == "pressure-dependent-Arrhenius": - for pressure, rate in reac.rate.rates: - reaction_lines.append( - PLOG_line.format( - pressure=pressure / ct.one_atm, - A=rate.pre_exponential_factor * unit_conversion_factor, - b=rate.temperature_exponent, - E_a=rate.activation_energy / CALORIES_CONSTANT, - ) - ) - elif reac.reaction_type == "chemically-activated": - rate = reac.rate + + elif reac.rate.type == "pressure-dependent-Arrhenius": + rate = reac.rate.rates[-1][1] reaction_lines.append( arrhenius_line.format( equation=reac.equation, max_reaction_length=max_reaction_length, - A=rate.low_rate.pre_exponential_factor * unit_conversion_factor, - b=rate.low_rate.temperature_exponent, - E_a=rate.low_rate.activation_energy / CALORIES_CONSTANT, - ) - ) - unit_conversion_factor /= 1_000.0 - reaction_lines.append( - high_line.format( - A=rate.high_rate.pre_exponential_factor * unit_conversion_factor, - b=rate.high_rate.temperature_exponent, - E_a=rate.high_rate.activation_energy / CALORIES_CONSTANT, + A=rate.pre_exponential_factor * unit_conversion_factor, + b=rate.temperature_exponent, + E_a=rate.activation_energy / CALORIES_CONSTANT, ) ) - if rate.type == "Troe": - reaction_lines.append( - "TROE /" - + " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs)) - + "/" - ) - elif rate.type == "SRI": + for pressure, rate in reac.rate.rates: reaction_lines.append( - "SRI /" - + " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs)) - + "/" + PLOG_line.format( + pressure=pressure / ct.one_atm, + A=rate.pre_exponential_factor * unit_conversion_factor, + b=rate.temperature_exponent, + E_a=rate.activation_energy / CALORIES_CONSTANT, + ) ) - elif reac.reaction_type == "falloff": + + elif reac.rate.type in ["falloff", "chemically-activated"]: rate = reac.rate + if reac.rate.type == "falloff": + rate1 = rate.high_rate + rate2 = rate.low_rate + unit_conversion_factor2 = unit_conversion_factor * 1_000.0 + other_line = low_line + else: + rate1 = rate.low_rate + rate2 = rate.high_rate + unit_conversion_factor2 = unit_conversion_factor / 1_000.0 + other_line = high_line + reaction_lines.append( arrhenius_line.format( equation=reac.equation, max_reaction_length=max_reaction_length, - A=rate.high_rate.pre_exponential_factor * unit_conversion_factor, - b=rate.high_rate.temperature_exponent, - E_a=rate.high_rate.activation_energy / CALORIES_CONSTANT, + A=rate1.pre_exponential_factor * unit_conversion_factor, + b=rate1.temperature_exponent, + E_a=rate1.activation_energy / CALORIES_CONSTANT, ) ) - unit_conversion_factor *= 1_000.0 reaction_lines.append( - low_line.format( - A=rate.low_rate.pre_exponential_factor * unit_conversion_factor, - b=rate.low_rate.temperature_exponent, - E_a=rate.low_rate.activation_energy / CALORIES_CONSTANT, + other_line.format( + A=rate2.pre_exponential_factor * unit_conversion_factor2, + b=rate2.temperature_exponent, + E_a=rate2.activation_energy / CALORIES_CONSTANT, ) ) - if rate.type == "Troe": + + if reac.rate.sub_type == "Troe": reaction_lines.append( "TROE /" + " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs)) + "/" ) - elif rate.type == "SRI": + elif reac.rate.sub_type == "SRI": reaction_lines.append( "SRI /" + " ".join(map(lambda s: format(s, ".7G"), rate.falloff_coeffs)) + "/" ) + elif reac.rate.sub_type != "Lindemann": + raise ValueError( + f"Unable to convert reaction type: '{reac.reaction_type}'" + ) + else: raise ValueError(f"Unknown reaction type: '{reac.reaction_type}'") - if getattr(reac, "efficiencies", None) is not None: + if reac.third_body is not None: reaction_lines.append( " ".join( - f"{spec}/{value:.3E}/" for spec, value in reac.efficiencies.items() + f"{spec}/{value:.3E}/" + for spec, value in reac.third_body.efficiencies.items() ) ) diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index 30e4e4bf47..8eda9b77cc 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -79,7 +79,9 @@ void ArrheniusBase::getRateParameters(AnyMap& node) const if (!valid()) { // Return empty/unmodified AnyMap return; - } else if (m_rate_units.factor() != 0.0) { + } + + if (m_rate_units.factor() != 0.0) { node[m_A_str].setQuantity(m_A, m_rate_units); } else { node[m_A_str] = m_A; @@ -117,9 +119,7 @@ void ArrheniusBase::getParameters(AnyMap& node) const { // RateType object is configured node["rate-constant"] = std::move(rateNode); } - if (type() != "Arrhenius") { - node["type"] = type(); - } + node["type"] = type(); } void ArrheniusBase::check(const std::string& equation) diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index 1a44e29b03..9ec80a5f30 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -136,9 +136,14 @@ bool BulkKinetics::addReaction(shared_ptr r, bool resize) } shared_ptr rate = r->rate(); + std::string rtype = rate->subType(); + if (rtype == "") { + rtype = rate->type(); + } + // If necessary, add new MultiRate evaluator - if (m_bulk_types.find(rate->type()) == m_bulk_types.end()) { - m_bulk_types[rate->type()] = m_bulk_rates.size(); + if (m_bulk_types.find(rtype) == m_bulk_types.end()) { + m_bulk_types[rtype] = m_bulk_rates.size(); m_bulk_rates.push_back(rate->newMultiRate()); m_bulk_rates.back()->resize(m_kk, nReactions(), nPhases()); } @@ -148,7 +153,7 @@ bool BulkKinetics::addReaction(shared_ptr r, bool resize) rate->setContext(*r, *this); // Add reaction rate to evaluator - size_t index = m_bulk_types[rate->type()]; + size_t index = m_bulk_types[rtype]; m_bulk_rates[index]->add(nReactions() - 1, *rate); // Add reaction to third-body evaluator @@ -186,14 +191,19 @@ void BulkKinetics::modifyReaction(size_t i, shared_ptr rNew) Kinetics::modifyReaction(i, rNew); shared_ptr rate = rNew->rate(); + std::string rtype = rate->subType(); + if (rtype == "") { + rtype = rate->type(); + } + // Ensure that MultiRate evaluator is available - if (m_bulk_types.find(rate->type()) == m_bulk_types.end()) { + if (m_bulk_types.find(rtype) == m_bulk_types.end()) { throw CanteraError("BulkKinetics::modifyReaction", - "Evaluator not available for type '{}'.", rate->type()); + "Evaluator not available for type '{}'.", rtype); } // Replace reaction rate to evaluator - size_t index = m_bulk_types[rate->type()]; + size_t index = m_bulk_types[rtype]; rate->setRateIndex(i); rate->setContext(*rNew, *this); diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 884847f341..0e0f9be695 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -428,15 +428,20 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base, bool resize) rate->setRateIndex(nReactions() - 1); rate->setContext(*r_base, *this); + std::string rtype = rate->subType(); + if (rtype == "") { + rtype = rate->type(); + } + // If necessary, add new interface MultiRate evaluator - if (m_interfaceTypes.find(rate->type()) == m_interfaceTypes.end()) { - m_interfaceTypes[rate->type()] = m_interfaceRates.size(); + if (m_interfaceTypes.find(rtype) == m_interfaceTypes.end()) { + m_interfaceTypes[rtype] = m_interfaceRates.size(); m_interfaceRates.push_back(rate->newMultiRate()); m_interfaceRates.back()->resize(m_kk, nReactions(), nPhases()); } // Add reaction rate to evaluator - size_t index = m_interfaceTypes[rate->type()]; + size_t index = m_interfaceTypes[rtype]; m_interfaceRates[index]->add(nReactions() - 1, *rate); return true; @@ -450,7 +455,11 @@ void InterfaceKinetics::modifyReaction(size_t i, shared_ptr r_base) rate->setRateIndex(i); rate->setContext(*r_base, *this); - const auto& rtype = rate->type(); + std::string rtype = rate->subType(); + if (rtype == "") { + rtype = rate->type(); + } + // Ensure that interface MultiRate evaluator is available if (!m_interfaceTypes.count(rtype)) { throw CanteraError("InterfaceKinetics::modifyReaction", diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index e17a5a6692..90b1cd7dc8 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -150,25 +150,9 @@ std::pair Kinetics::checkDuplicates(bool throw_err) const continue; // stoichiometries differ (not by a multiple) } else if (c < 0.0 && !R.reversible && !other.reversible) { continue; // irreversible reactions in opposite directions - } else if (R.type() == "falloff" || R.type() == "chemically-activated") { - auto tb1 = dynamic_cast(R).thirdBody(); - auto tb2 = dynamic_cast(other).thirdBody(); - bool thirdBodyOk = true; - for (size_t k = 0; k < nTotalSpecies(); k++) { - string s = kineticsSpeciesName(k); - if (tb1->efficiency(s) * tb2->efficiency(s) != 0.0) { - // non-zero third body efficiencies for species `s` in - // both reactions - thirdBodyOk = false; - break; - } - } - if (thirdBodyOk) { - continue; // No overlap in third body efficiencies - } - } else if (R.type() == "three-body") { - ThirdBody& tb1 = *(dynamic_cast(R).thirdBody()); - ThirdBody& tb2 = *(dynamic_cast(other).thirdBody()); + } else if (R.usesThirdBody() && other.usesThirdBody()) { + ThirdBody& tb1 = *(R.thirdBody()); + ThirdBody& tb2 = *(other.thirdBody()); bool thirdBodyOk = true; for (size_t k = 0; k < nTotalSpecies(); k++) { string s = kineticsSpeciesName(k); diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 77c0faaa30..87d5252ca3 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -6,9 +6,10 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/kinetics/Reaction.h" -#include "ReactionFactory.h" #include "cantera/kinetics/ReactionRateFactory.h" #include "cantera/kinetics/Kinetics.h" +#include "cantera/kinetics/Arrhenius.h" // @todo: remove after Cantera 3.0 +#include "cantera/kinetics/Falloff.h" // @todo: remove after Cantera 3.0 #include "cantera/thermo/ThermoPhase.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/base/Array.h" @@ -28,45 +29,52 @@ namespace Cantera Reaction::Reaction(const Composition& reactants_, const Composition& products_, - shared_ptr rate_) + shared_ptr rate_, + shared_ptr tbody_) : reactants(reactants_) , products(products_) - , m_rate(rate_) + , m_from_composition(true) + , m_third_body(tbody_) { + setRate(rate_); +} + +Reaction::Reaction(const std::string& equation, + shared_ptr rate_, + shared_ptr tbody_) + : m_third_body(tbody_) +{ + setEquation(equation); + setRate(rate_); } Reaction::Reaction(const AnyMap& node, const Kinetics& kin) - : Reaction() { + std::string rate_type = node.getString("type", "Arrhenius"); + if (!kin.nPhases()) { + throw InputFileError("Reaction", node, + "Cannot instantiate Reaction with empty Kinetics object."); + } + setParameters(node, kin); - if (kin.nPhases()) { - size_t nDim = kin.thermo(kin.reactionPhaseIndex()).nDim(); - if (nDim == 3) { - setRate(newReactionRate(node, calculateRateCoeffUnits(kin))); - } else { - AnyMap rateNode = node; - if (!rateNode.hasKey("type")) { - // Reaction type is not specified - rateNode["type"] = "Arrhenius"; + size_t nDim = kin.thermo(kin.reactionPhaseIndex()).nDim(); + if (nDim == 3) { + setRate(newReactionRate(node, calculateRateCoeffUnits(kin))); + } else { + AnyMap rateNode = node; + if (rateNode.hasKey("rate-constant")) { + if (!ba::starts_with(rate_type, "interface-")) { + rateNode["type"] = "interface-" + rate_type; } - std::string type = rateNode["type"].asString(); - if (rateNode.hasKey("rate-constant")) { - if (!boost::algorithm::starts_with(type, "interface-")) { - rateNode["type"] = "interface-" + type; - } - } else if (node.hasKey("sticking-coefficient")) { - if (!boost::algorithm::starts_with(type, "sticking-")) { - rateNode["type"] = "sticking-" + type; - } - } else { - throw InputFileError("Reaction::Reaction", input, - "Unable to infer interface reaction type."); + } else if (node.hasKey("sticking-coefficient")) { + if (!ba::starts_with(rate_type, "sticking-")) { + rateNode["type"] = "sticking-" + rate_type; } - setRate(newReactionRate(rateNode, calculateRateCoeffUnits(kin))); + } else { + throw InputFileError("Reaction::Reaction", input, + "Unable to infer interface reaction type."); } - } else { - // This route is used by the Python API. - setRate(newReactionRate(node)); + setRate(newReactionRate(rateNode, calculateRateCoeffUnits(kin))); } check(); } @@ -131,6 +139,11 @@ AnyMap Reaction::parameters(bool withInput) const void Reaction::getParameters(AnyMap& reactionNode) const { + if (!m_rate) { + throw CanteraError("Reaction::getParameters", + "Serialization of empty Reaction object is not supported."); + } + reactionNode["equation"] = equation(); if (duplicate) { @@ -146,26 +159,35 @@ void Reaction::getParameters(AnyMap& reactionNode) const reactionNode["nonreactant-orders"] = true; } - if (m_rate) { - reactionNode.update(m_rate->parameters()); + reactionNode.update(m_rate->parameters()); - // strip information not needed for reconstruction - if (reactionNode.hasKey("type")) { - std::string type = reactionNode["type"].asString(); - if (boost::algorithm::starts_with(type, "Arrhenius")) { - reactionNode.erase("type"); - } else if (boost::algorithm::starts_with(type, "Blowers-Masel")) { - reactionNode["type"] = "Blowers-Masel"; - } + // strip information not needed for reconstruction + std::string type = reactionNode["type"].asString(); + if (type == "pressure-dependent-Arrhenius") { + // skip + } else if (m_explicit_rate && ba::ends_with(type, "Arrhenius")) { + // retain type information + if (m_third_body) { + reactionNode["type"] = "three-body"; + } else { + reactionNode["type"] = "elementary"; } + } else if (ba::ends_with(type, "Arrhenius")) { + reactionNode.erase("type"); + } else if (ba::ends_with(type, "Blowers-Masel")) { + reactionNode["type"] = "Blowers-Masel"; + } + + if (m_third_body) { + m_third_body->getParameters(reactionNode); } } void Reaction::setParameters(const AnyMap& node, const Kinetics& kin) { if (node.empty()) { - // empty node: used by newReaction() factory loader - return; + throw InputFileError("Reaction::setParameters", input, + "Cannot set reaction parameters from empty node."); } input = node; @@ -186,28 +208,55 @@ void Reaction::setParameters(const AnyMap& node, const Kinetics& kin) duplicate = node.getBool("duplicate", false); allow_negative_orders = node.getBool("negative-orders", false); allow_nonreactant_orders = node.getBool("nonreactant-orders", false); + + if (m_third_body) { + m_third_body->setParameters(node); + } else if (node.hasKey("default-efficiency") || node.hasKey("efficiencies")) { + throw InputFileError("Reaction::setParameters", input, + "Reaction '{}' specifies efficiency parameters\n" + "but does not involve third body colliders.", equation()); + } } void Reaction::setRate(shared_ptr rate) { if (!rate) { - // null pointer - m_rate.reset(); - } else { - m_rate = rate; - } - - if (reactants.count("(+M)") && std::dynamic_pointer_cast(m_rate)) { - warn_deprecated("Chebyshev reaction equation", input, "Specifying '(+M)' " - "in the reaction equation for Chebyshev reactions is deprecated."); - // remove optional third body notation - reactants.erase("(+M)"); - products.erase("(+M)"); + throw InputFileError("Reaction::setRate", input, + "Reaction rate for reaction '{}' must not be empty.", equation()); } + m_rate = rate; - if (reactants.count("M") && std::dynamic_pointer_cast(m_rate)) { - throw InputFileError("Reaction::setRate", input, "Found superfluous 'M' in " - "pressure-dependent-Arrhenius reaction."); + std::string rate_type = m_rate->type(); + if (m_third_body) { + if (rate_type == "falloff" || rate_type == "chemically-activated") { + if (m_third_body->mass_action && !m_from_composition) { + throw InputFileError("Reaction::setRate", input, + "Third-body collider does not use '(+{})' notation.", + m_third_body->name()); + } + m_third_body->mass_action = false; + } else if (rate_type == "Chebyshev") { + if (m_third_body->name() == "M") { + warn_deprecated("Chebyshev reaction equation", input, "Specifying 'M' " + "in the reaction equation for Chebyshev reactions is deprecated."); + m_third_body.reset(); + } + } else if (rate_type == "pressure-dependent-Arrhenius") { + if (m_third_body->name() == "M") { + throw InputFileError("Reaction::setRate", input, + "Found superfluous '{}' in pressure-dependent-Arrhenius reaction.", + m_third_body->name()); + } + } + } else { + if (rate_type == "falloff" || rate_type == "chemically-activated") { + if (!m_from_composition) { + throw InputFileError("Reaction::setRate", input, + "Reaction equation for falloff reaction '{}'\n does not " + "contain valid pressure-dependent third body", equation()); + } + m_third_body.reset(new ThirdBody("(+M)")); + } } } @@ -223,7 +272,10 @@ std::string Reaction::reactantString() const } result << iter->first; } - return result.str(); + if (m_third_body) { + result << m_third_body->collider(); + } + return result.str(); } std::string Reaction::productString() const @@ -238,7 +290,10 @@ std::string Reaction::productString() const } result << iter->first; } - return result.str(); + if (m_third_body) { + result << m_third_body->collider(); + } + return result.str(); } std::string Reaction::equation() const @@ -253,11 +308,127 @@ std::string Reaction::equation() const void Reaction::setEquation(const std::string& equation, const Kinetics* kin) { parseReactionEquation(*this, equation, input, kin); + + std::string rate_type = input.getString("type", ""); + if (rate_type == "three-body") { + // state type when serializing + m_explicit_rate = true; + } else if (rate_type == "elementary") { + // user override + m_explicit_rate = true; + return; + } else if (kin && kin->thermo(kin->reactionPhaseIndex()).nDim() != 3) { + // interface reactions + return; + } + + std::string third_body; + size_t count = 0; + size_t countM = 0; + for (const auto& reac : reactants) { + // detect explicitly specified collision partner + if (products.count(reac.first)) { + third_body = reac.first; + size_t generic = third_body == "M" + || third_body == "(+M)" || third_body == "(+ M)"; + count++; + countM += generic; + if (reac.second > 1 && products[third_body] > 1) { + count++; + countM += generic; + } + } + } + + if (count == 0) { + if (rate_type == "three-body") { + throw InputFileError("Reaction::setEquation", input, + "Reactants for reaction '{}'\n" + "do not contain a valid third body collider", equation); + } + return; + + } else if (countM > 1) { + throw InputFileError("Reaction::setEquation", input, + "Multiple generic third body colliders 'M' are not supported", equation); + + } else if (count > 1) { + // equations with more than one explicit third-body collider are handled as a + // regular elementary reaction unless the equation contains a generic third body + if (!countM) { + return; + } + third_body = "M"; + + } else if (!ba::starts_with(third_body, "(+")) { + // check for conditions of three-body reactions: + // - integer stoichiometric conditions + // - either reactant or product side involves exactly three species + size_t nreac = 0; + size_t nprod = 0; + + // ensure that all reactants have integer stoichiometric coefficients + for (const auto& reac : reactants) { + if (trunc(reac.second) != reac.second) { + return; + } + nreac += static_cast(reac.second); + } + + // ensure that all products have integer stoichiometric coefficients + for (const auto& prod : products) { + if (trunc(prod.second) != prod.second) { + return; + } + nprod += static_cast(prod.second); + } + + // either reactant or product side involves exactly three species + if (nreac != 3 && nprod != 3) { + return; + } + } + + if (m_third_body) { + m_third_body->setName(third_body); + } else { + m_third_body.reset(new ThirdBody(third_body)); + } + + // adjust reactant coefficients + auto reac = reactants.find(third_body); + if (trunc(reac->second) != 1) { + reac->second -= 1.; + } else { + reactants.erase(reac); + } + + // adjust product coefficients + auto prod = products.find(third_body); + if (trunc(prod->second) != 1) { + prod->second -= 1.; + } else { + products.erase(prod); + } } std::string Reaction::type() const { - return "reaction"; + if (!m_rate) { + throw CanteraError("Reaction::type", "Empty Reaction does not have a type"); + } + + std::string rate_type = m_rate->type(); + std::string sub_type = m_rate->subType(); + if (sub_type != "") { + return rate_type + "-" + sub_type; + } + + if (m_third_body) { + return "three-body-" + rate_type; + } + + return rate_type; } UnitStack Reaction::calculateRateCoeffUnits(const Kinetics& kin) @@ -313,17 +484,6 @@ void updateUndeclared(std::vector& undeclared, } } -std::pair, bool> Reaction::undeclaredThirdBodies( - const Kinetics& kin) const -{ - std::vector undeclared; - if (m_third_body) { - updateUndeclared(undeclared, m_third_body->efficiencies, kin); - return std::make_pair(undeclared, m_third_body->specified_collision_partner); - } - return std::make_pair(undeclared, false); -} - void Reaction::checkBalance(const Kinetics& kin) const { Composition balr, balp; @@ -406,7 +566,7 @@ bool Reaction::checkSpecies(const Kinetics& kin) const } else { throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n" "contains undeclared species: '{}'", - equation(), boost::algorithm::join(undeclared, "', '")); + equation(), ba::join(undeclared, "', '")); } } @@ -420,34 +580,17 @@ bool Reaction::checkSpecies(const Kinetics& kin) const throw InputFileError("Reaction::checkSpecies", input["orders"], "Reaction '{}'\n" "defines reaction orders for undeclared species: '{}'", - equation(), boost::algorithm::join(undeclared, "', '")); + equation(), ba::join(undeclared, "', '")); } // Error for empty input AnyMap (that is, XML) throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n" "defines reaction orders for undeclared species: '{}'", - equation(), boost::algorithm::join(undeclared, "', '")); + equation(), ba::join(undeclared, "', '")); } } - // Use helper function while there is no uniform handling of third bodies - auto third = undeclaredThirdBodies(kin); - undeclared = third.first; - bool specified_collision_partner_ = third.second; - if (!undeclared.empty()) { - if (!kin.skipUndeclaredThirdBodies()) { - if (input.hasKey("efficiencies")) { - throw InputFileError("Reaction::checkSpecies", input["efficiencies"], - "Reaction '{}'\n" - "defines third-body efficiencies for undeclared species: '{}'", - equation(), boost::algorithm::join(undeclared, "', '")); - } - // Error for specified ThirdBody or empty input AnyMap - throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n" - "is a three-body reaction with undeclared species: '{}'", - equation(), boost::algorithm::join(undeclared, "', '")); - } else if (kin.skipUndeclaredSpecies() && specified_collision_partner_) { - return false; - } + if (m_third_body) { + return m_third_body->checkSpecies(*this, kin); } checkBalance(kin); @@ -488,365 +631,188 @@ bool Reaction::usesElectrochemistry(const Kinetics& kin) const ThirdBody::ThirdBody(double default_eff) : default_efficiency(default_eff) - , specified_collision_partner(false) - , mass_action(true) { + warn_deprecated("ThirdBody", + "Instantiation with default efficiency is deprecated and will be removed " + "after Cantera 3.0. Instantiate with collider name instead."); } -ThirdBody::ThirdBody(const AnyMap& node) - : specified_collision_partner(false) - , mass_action(true) +ThirdBody::ThirdBody(const std::string& third_body) { - setEfficiencies(node); + setName(third_body); } -void ThirdBody::setEfficiencies(const AnyMap& node) +void ThirdBody::setName(const std::string& third_body) { - default_efficiency = node.getDouble("default-efficiency", 1.0); - if (node.hasKey("efficiencies")) { - efficiencies = node["efficiencies"].asMap(); + std::string name = third_body; + if (ba::starts_with(third_body, "(+ ")) { + mass_action = false; + name = third_body.substr(3, third_body.size() - 4); + } else if (ba::starts_with(third_body, "(+")) { + mass_action = false; + name = third_body.substr(2, third_body.size() - 3); } -} -double ThirdBody::efficiency(const std::string& k) const -{ - return getValue(efficiencies, k, default_efficiency); + if (name == m_name) { + return; + } + if (name == "M") { + throw CanteraError("ThirdBody::setName", + "Unable to revert explicit third body '{}' to 'M'", m_name); + } + if (efficiencies.size()) { + throw CanteraError("ThirdBody::setName", + "Conflicting efficiency definition for explicit third body '{}'", name); + } + m_name = name; + default_efficiency = 0.; + efficiencies[m_name] = 1.; } -ThreeBodyReaction::ThreeBodyReaction() +ThirdBody::ThirdBody(const AnyMap& node) { - m_third_body.reset(new ThirdBody); - setRate(newReactionRate(type())); + setParameters(node); } -ThreeBodyReaction::ThreeBodyReaction(const Composition& reactants, - const Composition& products, - const ArrheniusRate& rate, - const ThirdBody& tbody) - : Reaction(reactants, products, make_shared(rate)) +void ThirdBody::setEfficiencies(const AnyMap& node) { - m_third_body = std::make_shared(tbody); + warn_deprecated("ThirdBody::setEfficiencies", node, + "To be removed after Cantera 3.0. Renamed to setParameters"); + setParameters(node); } -ThreeBodyReaction::ThreeBodyReaction(const AnyMap& node, const Kinetics& kin) +void ThirdBody::setParameters(const AnyMap& node) { - m_third_body.reset(new ThirdBody); - if (!node.empty()) { - setParameters(node, kin); - setRate(newReactionRate(node, calculateRateCoeffUnits(kin))); - } else { - setRate(newReactionRate(type())); + if (node.hasKey("default-efficiency")) { + default_efficiency = node["default-efficiency"].asDouble(); + } + if (node.hasKey("efficiencies")) { + efficiencies = node["efficiencies"].asMap(); } } -bool ThreeBodyReaction::detectEfficiencies() +void ThirdBody::getParameters(AnyMap& node) const { - for (const auto& reac : reactants) { - // detect explicitly specified collision partner - if (products.count(reac.first)) { - m_third_body->efficiencies[reac.first] = 1.; + if (m_name == "M") { + if (efficiencies.size()) { + node["efficiencies"] = efficiencies; + node["efficiencies"].setFlowStyle(); + } + if (default_efficiency != 1.0) { + node["default-efficiency"] = default_efficiency; } } - - if (m_third_body->efficiencies.size() == 0) { - return false; - } else if (m_third_body->efficiencies.size() > 1) { - throw CanteraError("ThreeBodyReaction::detectEfficiencies", - "Found more than one explicitly specified collision partner\n" - "in reaction '{}'.", equation()); - } - - m_third_body->default_efficiency = 0.; - m_third_body->specified_collision_partner = true; - auto sp = m_third_body->efficiencies.begin(); - - // adjust reactant coefficients - auto reac = reactants.find(sp->first); - if (trunc(reac->second) != 1) { - reac->second -= 1.; - } else { - reactants.erase(reac); - } - - // adjust product coefficients - auto prod = products.find(sp->first); - if (trunc(prod->second) != 1) { - prod->second -= 1.; - } else { - products.erase(prod); - } - - return true; } -void ThreeBodyReaction::setEquation(const std::string& equation, const Kinetics* kin) +double ThirdBody::efficiency(const std::string& k) const { - Reaction::setEquation(equation, kin); - if (reactants.count("M") != 1 || products.count("M") != 1) { - if (!detectEfficiencies()) { - throw InputFileError("ThreeBodyReaction::setParameters", input, - "Reaction equation '{}' does not contain third body 'M'", - equation); - } - return; - } - - reactants.erase("M"); - products.erase("M"); + return getValue(efficiencies, k, default_efficiency); } -void ThreeBodyReaction::setParameters(const AnyMap& node, const Kinetics& kin) +std::string ThirdBody::collider() const { - if (node.empty()) { - // empty node: used by newReaction() factory loader - return; - } - Reaction::setParameters(node, kin); - if (!m_third_body->specified_collision_partner) { - m_third_body->setEfficiencies(node); + if (mass_action) { + return " + " + m_name; } + return " (+" + m_name + ")"; } -void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const +bool ThirdBody::checkSpecies(const Reaction& rxn, const Kinetics& kin) const { - Reaction::getParameters(reactionNode); - if (!m_third_body->specified_collision_partner) { - reactionNode["type"] = "three-body"; - reactionNode["efficiencies"] = m_third_body->efficiencies; - reactionNode["efficiencies"].setFlowStyle(); - if (m_third_body->default_efficiency != 1.0) { - reactionNode["default-efficiency"] = m_third_body->default_efficiency; + std::vector undeclared; + updateUndeclared(undeclared, efficiencies, kin); + + if (!undeclared.empty()) { + if (!kin.skipUndeclaredThirdBodies()) { + if (rxn.input.hasKey("efficiencies")) { + throw InputFileError("ThirdBody::checkSpecies", + rxn.input["efficiencies"], "Reaction '{}'\n" + "defines third-body efficiencies for undeclared species: '{}'", + rxn.equation(), ba::join(undeclared, "', '")); + } + // Error for specified ThirdBody or empty input AnyMap + throw InputFileError("ThirdBody::checkSpecies", rxn.input, "Reaction '{}'\n" + "is a three-body reaction with undeclared species: '{}'", + rxn.equation(), ba::join(undeclared, "', '")); + } else if (kin.skipUndeclaredSpecies() && m_name != "M") { + return false; } } + return true; } -std::string ThreeBodyReaction::reactantString() const +ThreeBodyReaction::ThreeBodyReaction() { - if (m_third_body->specified_collision_partner) { - return Reaction::reactantString() + " + " - + m_third_body->efficiencies.begin()->first; - } else { - return Reaction::reactantString() + " + M"; - } + warn_deprecated("ThreeBodyReaction", + "To be removed after Cantera 3.0. Replaceable with Reaction."); + m_third_body.reset(new ThirdBody); + setRate(newReactionRate(type())); } -std::string ThreeBodyReaction::productString() const +ThreeBodyReaction::ThreeBodyReaction(const Composition& reactants, + const Composition& products, + const ArrheniusRate& rate, + const ThirdBody& tbody) + : Reaction(reactants, + products, + make_shared(rate), + make_shared(tbody)) { - if (m_third_body->specified_collision_partner) { - return Reaction::productString() + " + " - + m_third_body->efficiencies.begin()->first; - } else { - return Reaction::productString() + " + M"; - } + warn_deprecated("ThreeBodyReaction", + "To be removed after Cantera 3.0. Replaceable with Reaction."); +} + +ThreeBodyReaction::ThreeBodyReaction(const AnyMap& node, const Kinetics& kin) + : Reaction(node, kin) +{ + warn_deprecated("ThreeBodyReaction", + "To be removed after Cantera 3.0. Replaceable with Reaction."); } FalloffReaction::FalloffReaction() { + warn_deprecated("FalloffReaction", + "To be removed after Cantera 3.0. Replaceable with Reaction."); m_third_body.reset(new ThirdBody); - m_third_body->mass_action = false; setRate(newReactionRate(type())); } -FalloffReaction::FalloffReaction(const Composition& reactants, - const Composition& products, - const ReactionRate& rate, - const ThirdBody& tbody) - : Reaction(reactants, products) +FalloffReaction::FalloffReaction(const Composition& reactants_, + const Composition& products_, + const FalloffRate& rate_, + const ThirdBody& tbody_) { - m_third_body = std::make_shared(tbody); - m_third_body->mass_action = false; - AnyMap node = rate.parameters(); + warn_deprecated("FalloffReaction", + "To be removed after Cantera 3.0. Replaceable with Reaction."); + // cannot be delegated as std::make_shared does not work for FalloffRate + reactants = reactants_; + products = products_; + m_third_body = std::make_shared(tbody_); + AnyMap node = rate_.parameters(); node.applyUnits(); std::string rate_type = node["type"].asString(); - if (rate_type != "falloff" && rate_type != "chemically-activated") { - // use node information to determine whether rate is a falloff rate - throw CanteraError("FalloffReaction::FalloffReaction", - "Incompatible types: '{}' is not a falloff rate object.", rate.type()); - } setRate(newReactionRate(node)); } FalloffReaction::FalloffReaction(const AnyMap& node, const Kinetics& kin) + : Reaction(node, kin) { - m_third_body.reset(new ThirdBody); - m_third_body->mass_action = false; - if (!node.empty()) { - setParameters(node, kin); - setRate(newReactionRate(node, calculateRateCoeffUnits(kin))); - } else { - setRate(newReactionRate(type())); - } -} - -std::string FalloffReaction::type() const -{ - if (m_rate && - std::dynamic_pointer_cast(m_rate)->chemicallyActivated()) - { - return "chemically-activated"; - } - return "falloff"; -} - -std::string FalloffReaction::reactantString() const -{ - if (m_third_body->specified_collision_partner) { - return Reaction::reactantString() + " (+" + - m_third_body->efficiencies.begin()->first + ")"; - } else { - return Reaction::reactantString() + " (+M)"; - } -} - -std::string FalloffReaction::productString() const -{ - if (m_third_body->specified_collision_partner) { - return Reaction::productString() + " (+" + - m_third_body->efficiencies.begin()->first + ")"; - } else { - return Reaction::productString() + " (+M)"; - } -} - -void FalloffReaction::setParameters(const AnyMap& node, const Kinetics& kin) -{ + warn_deprecated("FalloffReaction", + "To be removed after Cantera 3.0. Replaceable with Reaction."); if (node.empty()) { - // empty node: used by newReaction() factory loader - return; - } - Reaction::setParameters(node, kin); - - if (!m_third_body->specified_collision_partner) { - m_third_body->setEfficiencies(node); - } -} - -void FalloffReaction::setEquation(const std::string& equation, const Kinetics* kin) -{ - Reaction::setEquation(equation, kin); - - // Detect falloff third body based on partial setup; - // parseReactionEquation (called via Reaction::setEquation) sets the - // stoichiometric coefficient of the falloff species to -1. - std::string third_body_str; - std::string third_body; - for (auto& reactant : reactants) { - if (reactant.second == -1 && ba::starts_with(reactant.first, "(+")) { - third_body_str = reactant.first; - third_body = third_body_str.substr(2, third_body_str.size() - 3); - break; - } - } - - // Equation must contain a third body, and it must appear on both sides - if (third_body_str == "") { - throw InputFileError("FalloffReaction::setParameters", input, - "Reactants for reaction '{}' do not contain a pressure-dependent " - "third body", equation); - } - if (products.count(third_body_str) == 0) { - throw InputFileError("FalloffReaction::setParameters", input, - "Unable to match third body '{}' in reactants and products of " - "reaction '{}'", third_body, equation); - } - - // Remove the dummy species - reactants.erase(third_body_str); - products.erase(third_body_str); - - if (third_body == "M") { - m_third_body->specified_collision_partner = false; - } else { - // Specific species is listed as the third body - m_third_body->default_efficiency = 0; - m_third_body->efficiencies.emplace(third_body, 1.0); - m_third_body->specified_collision_partner = true; - } -} - -void FalloffReaction::getParameters(AnyMap& reactionNode) const -{ - Reaction::getParameters(reactionNode); - if (m_third_body->specified_collision_partner) { - // pass - } else if (m_third_body->efficiencies.size()) { - reactionNode["efficiencies"] = m_third_body->efficiencies; - reactionNode["efficiencies"].setFlowStyle(); - if (m_third_body->default_efficiency != 1.0) { - reactionNode["default-efficiency"] = m_third_body->default_efficiency; - } - } -} - -bool isThreeBody(const Reaction& R) -{ - // detect explicitly specified collision partner - size_t found = 0; - for (const auto& reac : R.reactants) { - auto prod = R.products.find(reac.first); - if (prod != R.products.end() && - trunc(reac.second) == reac.second && trunc(prod->second) == prod->second) { - // candidate species with integer stoichiometric coefficients on both sides - found += 1; - } - } - if (found != 1) { - return false; - } - - // ensure that all reactants have integer stoichiometric coefficients - size_t nreac = 0; - for (const auto& reac : R.reactants) { - if (trunc(reac.second) != reac.second) { - return false; - } - nreac += static_cast(reac.second); - } - - // ensure that all products have integer stoichiometric coefficients - size_t nprod = 0; - for (const auto& prod : R.products) { - if (trunc(prod.second) != prod.second) { - return false; - } - nprod += static_cast(prod.second); + m_third_body.reset(new ThirdBody); + setRate(newReactionRate(type())); } - - // either reactant or product side involves exactly three species - return (nreac == 3) || (nprod == 3); } unique_ptr newReaction(const std::string& type) { - AnyMap rxn_node; - Kinetics kin; - unique_ptr R(ReactionFactory::factory()->create(type, rxn_node, kin)); - return R; + return unique_ptr(new Reaction()); } unique_ptr newReaction(const AnyMap& rxn_node, const Kinetics& kin) { - std::string type = "elementary"; - size_t nDim = kin.thermo(kin.reactionPhaseIndex()).nDim(); - if (rxn_node.hasKey("type")) { - type = rxn_node["type"].asString(); - } else if (nDim == 3) { - // Reaction type is not specified - // See if this is a three-body reaction with a specified collision partner - Reaction testReaction; - parseReactionEquation(testReaction, rxn_node["equation"].asString(), - rxn_node, &kin); - if (isThreeBody(testReaction)) { - type = "three-body"; - } - } - - if (!(ReactionFactory::factory()->exists(type))) { - throw InputFileError("ReactionFactory::newReaction", rxn_node["type"], - "Unknown reaction type '{}'", type); - } - Reaction* R = ReactionFactory::factory()->create(type, rxn_node, kin); - return unique_ptr(R); + return unique_ptr(new Reaction(rxn_node, kin)); } void parseReactionEquation(Reaction& R, const std::string& equation, @@ -865,18 +831,21 @@ void parseReactionEquation(Reaction& R, const std::string& equation, tokens[i] == "<=>" || tokens[i] == "=" || tokens[i] == "=>") { std::string species = tokens[i-1]; - double stoich; - if (last_used != npos && tokens[last_used] == "(+") { + double stoich = 1.0; + bool mass_action = true; + if (last_used != npos && tokens[last_used] == "(+" + && ba::ends_with(species, ")")) { // Falloff third body with space, such as "(+ M)" + mass_action = false; species = "(+" + species; - stoich = -1; - } else if (last_used == i-1 && ba::starts_with(species, "(+") - && ba::ends_with(species, ")")) { + } else if (last_used == i - 1 && ba::starts_with(species, "(+") + && ba::ends_with(species, ")")) { // Falloff 3rd body written without space, such as "(+M)" - stoich = -1; - } else if (last_used == i-2) { // Species with no stoich. coefficient - stoich = 1.0; - } else if (last_used == i-3) { // Stoich. coefficient and species + mass_action = false; + } else if (last_used == i - 2) { + // Species with no stoich. coefficient + } else if (last_used == i - 3) { + // Stoich. coefficient and species try { stoich = fpValueCheck(tokens[i-2]); } catch (CanteraError& err) { @@ -891,7 +860,7 @@ void parseReactionEquation(Reaction& R, const std::string& equation, (last_used == npos) ? "n/a" : tokens[last_used]); } if (!kin || (kin->kineticsSpeciesIndex(species) == npos - && stoich != -1 && species != "M")) + && mass_action && species != "M")) { R.setValid(false); } @@ -921,7 +890,7 @@ std::vector> getReactions(const AnyValue& items, { std::vector> all_reactions; for (const auto& node : items.asVector()) { - shared_ptr R(newReaction(node, kinetics)); + shared_ptr R(new Reaction(node, kinetics)); R->check(); R->validate(kinetics); if (R->valid() && R->checkSpecies(kinetics)) { diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp deleted file mode 100644 index 2b0bdd2c1f..0000000000 --- a/src/kinetics/ReactionFactory.cpp +++ /dev/null @@ -1,75 +0,0 @@ - /** - * @file ReactionFactory.cpp - */ - -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#include "cantera/thermo/ThermoPhase.h" -#include "cantera/kinetics/Reaction.h" -#include "ReactionFactory.h" -#include "cantera/kinetics/ReactionRateFactory.h" -#include "cantera/kinetics/Kinetics.h" -#include "cantera/base/AnyMap.h" -#include "cantera/base/stringUtils.h" - -namespace Cantera -{ - -ReactionFactory* ReactionFactory::s_factory = 0; -std::mutex ReactionFactory::reaction_mutex; - -ReactionFactory::ReactionFactory() -{ - // register elementary reactions - reg("reaction", [](const AnyMap& node, const Kinetics& kin) { - return new Reaction(node, kin); - }); - addAlias("reaction", "elementary"); - addAlias("reaction", "arrhenius"); - addAlias("reaction", "Arrhenius"); - addAlias("reaction", ""); - - // register three-body reactions - reg("three-body", [](const AnyMap& node, const Kinetics& kin) { - return new ThreeBodyReaction(node, kin); - }); - addAlias("three-body", "threebody"); - addAlias("three-body", "three_body"); - - // register falloff reactions - reg("falloff", [](const AnyMap& node, const Kinetics& kin) { - return new FalloffReaction(node, kin); - }); - addAlias("falloff", "chemically-activated"); - addAlias("falloff", "chemact"); - addAlias("falloff", "chemically_activated"); - - addAlias("reaction", "pressure-dependent-Arrhenius"); - addAlias("reaction", "plog"); - addAlias("reaction", "pdep_arrhenius"); - - // register Chebyshev reactions - addAlias("reaction", "Chebyshev"); - addAlias("reaction", "chebyshev"); - - // register custom reactions specified by Func1 objects - addAlias("reaction", "custom-rate-function"); - - addAlias("reaction", "interface-Arrhenius"); - addAlias("reaction", "sticking-Arrhenius"); - - // register legacy interface reaction names - addAlias("reaction", "interface"); - addAlias("reaction", "surface"); - addAlias("reaction", "edge"); - addAlias("reaction", "electrochemical"); - - addAlias("reaction", "two-temperature-plasma"); - - addAlias("reaction", "Blowers-Masel"); - addAlias("reaction", "interface-Blowers-Masel"); - addAlias("reaction", "sticking-Blowers-Masel"); -} - -} diff --git a/src/kinetics/ReactionFactory.h b/src/kinetics/ReactionFactory.h deleted file mode 100644 index a014fbb63d..0000000000 --- a/src/kinetics/ReactionFactory.h +++ /dev/null @@ -1,64 +0,0 @@ -/** - * @file ReactionFactory.h - * Factory class for reaction functions. Used by classes - * that implement kinetics - * (see \ref reactionGroup and class \link Cantera::Reaction Reaction\endlink). - */ - -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#ifndef CT_NEWREACTION_H -#define CT_NEWREACTION_H - -#include "cantera/base/FactoryBase.h" -#include "cantera/kinetics/Reaction.h" - -namespace Cantera -{ - -/** - * Factory class to construct reaction function calculators. - * The reaction factory is accessed through static method factory: - * - * @code - * Reaction* f = ReactionFactory::factory()->newReaction(type, c) - * @endcode - * - * @ingroup reactionGroup - */ -class ReactionFactory : public Factory -{ -public: - /** - * Return a pointer to the factory. On the first call, a new instance is - * created. Since there is no need to instantiate more than one factory, - * on all subsequent calls, a pointer to the existing factory is returned. - */ - static ReactionFactory* factory() { - std::unique_lock lock(reaction_mutex); - if (!s_factory) { - s_factory = new ReactionFactory; - } - return s_factory; - } - - virtual void deleteFactory() { - std::unique_lock lock(reaction_mutex); - delete s_factory; - s_factory = 0; - } - -private: - //! Pointer to the single instance of the factory - static ReactionFactory* s_factory; - - //! default constructor, which is defined as private - ReactionFactory(); - - //! Mutex for use when calling the factory - static std::mutex reaction_mutex; -}; - -} -#endif diff --git a/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 356d8e003a..3cd354b368 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -14,8 +14,8 @@ #include "cantera/kinetics/Custom.h" #include "cantera/kinetics/Falloff.h" #include "cantera/kinetics/InterfaceRate.h" -#include "cantera/kinetics/TwoTempPlasmaRate.h" #include "cantera/kinetics/PlogRate.h" +#include "cantera/kinetics/TwoTempPlasmaRate.h" namespace Cantera { diff --git a/test/data/kineticsfromscratch.yaml b/test/data/kineticsfromscratch.yaml index f1196cc866..72a4585bad 100644 --- a/test/data/kineticsfromscratch.yaml +++ b/test/data/kineticsfromscratch.yaml @@ -98,6 +98,9 @@ reactions: - equation: O + OH => O + OH # Reaction 13 type: two-temperature-plasma rate-constant: {A: 17283, b: -3.1} +- equation: O + H2 + M <=> H2O + M # Reaction 14 + type: Blowers-Masel + rate-constant: {A: 38700, b: 2.7, Ea0: 2.619184e4, w: 4.184e9} surface: - units: {length: cm, quantity: mol, activation-energy: J/mol} diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index 80b8f04262..d860c47f9a 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -5,6 +5,13 @@ #include "cantera/thermo/Species.h" #include "cantera/kinetics/GasKinetics.h" #include "cantera/kinetics/InterfaceKinetics.h" +#include "cantera/kinetics/Arrhenius.h" +#include "cantera/kinetics/ChebyshevRate.h" +#include "cantera/kinetics/Custom.h" +#include "cantera/kinetics/Falloff.h" +#include "cantera/kinetics/InterfaceRate.h" +#include "cantera/kinetics/PlogRate.h" +#include "cantera/kinetics/TwoTempPlasmaRate.h" #include "cantera/base/Array.h" #include "cantera/base/stringUtils.h" @@ -49,7 +56,7 @@ class KineticsFromScratch : public testing::Test } }; -TEST_F(KineticsFromScratch, add_elementary_reaction) +TEST_F(KineticsFromScratch, add_elementary_reaction1) { // reaction 0: // equation: O + H2 <=> H + OH # Reaction 1 @@ -63,7 +70,20 @@ TEST_F(KineticsFromScratch, add_elementary_reaction) check_rates(0); } -TEST_F(KineticsFromScratch, add_three_body_reaction) +TEST_F(KineticsFromScratch, add_elementary_reaction2) +{ + // reaction 0: + // equation: O + H2 <=> H + OH # Reaction 1 + // rate-constant: {A: 38.7, b: 2.7, Ea: 6260.0} + std::string equation = "O + H2 <=> H + OH"; + auto rate = make_shared(3.87e1, 2.7, 2.619184e+07); + auto R = make_shared(equation, rate); + + kin.addReaction(R); + check_rates(0); +} + +TEST_F(KineticsFromScratch, add_three_body_reaction1) { // reaction 1: // equation: 2 O + M <=> O2 + M # Reaction 2 @@ -72,23 +92,79 @@ TEST_F(KineticsFromScratch, add_three_body_reaction) // efficiencies: {AR: 0.83, H2: 2.4, H2O: 15.4} Composition reac = parseCompString("O:2"); Composition prod = parseCompString("O2:1"); - ArrheniusRate rate(1.2e11, -1.0, 0.0); - ThirdBody tbody; - tbody.efficiencies = parseCompString("AR:0.83 H2:2.4 H2O:15.4"); - auto R = make_shared(reac, prod, rate, tbody); + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto tbody = make_shared(); + tbody->efficiencies = parseCompString("AR:0.83 H2:2.4 H2O:15.4"); + auto R = make_shared(reac, prod, rate, tbody); + + kin.addReaction(R); + check_rates(1); +} + +TEST_F(KineticsFromScratch, add_three_body_reaction2) +{ + // reaction 1: + // equation: 2 O + M <=> O2 + M # Reaction 2 + // type: three-body + // rate-constant: {A: 1.2e+11, b: -1.0, Ea: 0.0} + // efficiencies: {AR: 0.83, H2: 2.4, H2O: 15.4} + std::string equation = "2 O + M <=> O2 + M"; + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto tbody = make_shared(); + tbody->efficiencies = parseCompString("AR:0.83 H2:2.4 H2O:15.4"); + auto R = make_shared(equation, rate, tbody); kin.addReaction(R); check_rates(1); } +TEST_F(KineticsFromScratch, add_three_body_reaction3) +{ + std::string equation = "2 O + M <=> O2 + M"; + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto R = make_shared(equation, rate); + EXPECT_TRUE(R->usesThirdBody()); +} + +TEST_F(KineticsFromScratch, multiple_third_bodies1) +{ + std::string equation = "2 H + 2 O2 <=> H2 + 2 O2"; + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto R = make_shared(equation, rate); + EXPECT_FALSE(R->usesThirdBody()); +} + +TEST_F(KineticsFromScratch, multiple_third_bodies2) +{ + std::string equation = "2 H + 2 M <=> H2 + 2 M"; + auto rate = make_shared(1.2e11, -1.0, 0.0); + ASSERT_THROW(Reaction(equation, rate), CanteraError); +} + +TEST_F(KineticsFromScratch, multiple_third_bodies3) +{ + std::string equation = "2 H + O2 + M <=> H2 + O2 + M"; + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto R = make_shared(equation, rate); + EXPECT_TRUE(R->usesThirdBody()); +} + +TEST_F(KineticsFromScratch, add_two_temperature_plasma) +{ + std::string equation = "O + H => O + H"; + auto rate = make_shared(17283, -3.1, -5820000, 1081000); + auto R = make_shared(equation, rate); + EXPECT_FALSE(R->usesThirdBody()); +} + TEST_F(KineticsFromScratch, undefined_third_body) { Composition reac = parseCompString("O:2"); Composition prod = parseCompString("O2:1"); - ArrheniusRate rate(1.2e11, -1.0, 0.0); - ThirdBody tbody; - tbody.efficiencies = parseCompString("H2:0.1 CO2:0.83"); - auto R = make_shared(reac, prod, rate, tbody); + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto tbody = make_shared(); + tbody->efficiencies = parseCompString("H2:0.1 CO2:0.83"); + auto R = make_shared(reac, prod, rate, tbody); ASSERT_THROW(kin.addReaction(R), CanteraError); } @@ -97,17 +173,17 @@ TEST_F(KineticsFromScratch, skip_undefined_third_body) { Composition reac = parseCompString("O:2"); Composition prod = parseCompString("O2:1"); - ArrheniusRate rate(1.2e11, -1.0, 0.0); - ThirdBody tbody; - tbody.efficiencies = parseCompString("H2:0.1 CO2:0.83"); - auto R = make_shared(reac, prod, rate, tbody); + auto rate = make_shared(1.2e11, -1.0, 0.0); + auto tbody = make_shared(); + tbody->efficiencies = parseCompString("H2:0.1 CO2:0.83"); + auto R = make_shared(reac, prod, rate, tbody); kin.skipUndeclaredThirdBodies(true); kin.addReaction(R); ASSERT_EQ((size_t) 1, kin.nReactions()); } -TEST_F(KineticsFromScratch, add_falloff_reaction) +TEST_F(KineticsFromScratch, add_falloff_reaction1) { // reaction 2: // equation: 2 OH (+ M) <=> H2O2 (+ M) # Reaction 3 @@ -122,14 +198,48 @@ TEST_F(KineticsFromScratch, add_falloff_reaction) ArrheniusRate high_rate(7.4e10, -0.37, 0.0); ArrheniusRate low_rate(2.3e12, -0.9, -7112800.0); vector_fp falloff_params { 0.7346, 94.0, 1756.0, 5182.0 }; - TroeRate rate(low_rate, high_rate, falloff_params); - ThirdBody tbody; - tbody.efficiencies = parseCompString("AR:0.7 H2:2.0 H2O:6.0"); - auto R = make_shared(reac, prod, rate, tbody); + auto rate = make_shared(low_rate, high_rate, falloff_params); + auto tbody = make_shared(); + tbody->efficiencies = parseCompString("AR:0.7 H2:2.0 H2O:6.0"); + auto R = make_shared(reac, prod, rate, tbody); + + kin.addReaction(R); + check_rates(2); +} + +TEST_F(KineticsFromScratch, add_falloff_reaction2) +{ + // reaction 2: + // equation: 2 OH (+ M) <=> H2O2 (+ M) # Reaction 3 + // duplicate: true + // type: falloff + // low-P-rate-constant: {A: 2.3e+12, b: -0.9, Ea: -1700.0} + // high-P-rate-constant: {A: 7.4e+10, b: -0.37, Ea: 0.0} + // Troe: {A: 0.7346, T3: 94.0, T1: 1756.0, T2: 5182.0} + // efficiencies: {AR: 0.7, H2: 2.0, H2O: 6.0} + std::string equation = "2 OH (+ M) <=> H2O2 (+ M)"; + ArrheniusRate high_rate(7.4e10, -0.37, 0.0); + ArrheniusRate low_rate(2.3e12, -0.9, -7112800.0); + vector_fp falloff_params { 0.7346, 94.0, 1756.0, 5182.0 }; + auto rate = make_shared(low_rate, high_rate, falloff_params); + auto tbody = make_shared(); + tbody->efficiencies = parseCompString("AR:0.7 H2:2.0 H2O:6.0"); + auto R = make_shared(equation, rate, tbody); + kin.addReaction(R); check_rates(2); } +TEST_F(KineticsFromScratch, missing_third_body) +{ + std::string equation = "2 OH <=> H2O2"; + ArrheniusRate high_rate(7.4e10, -0.37, 0.0); + ArrheniusRate low_rate(2.3e12, -0.9, -7112800.0); + vector_fp falloff_params { 0.7346, 94.0, 1756.0, 5182.0 }; + auto rate = make_shared(low_rate, high_rate, falloff_params); + ASSERT_THROW(Reaction(equation, rate), CanteraError); +} + TEST_F(KineticsFromScratch, add_plog_reaction) { // reaction 3: diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 7b81482832..703297ea82 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -2,10 +2,16 @@ #include "cantera/base/Units.h" #include "cantera/base/Solution.h" #include "cantera/base/Interface.h" +#include "cantera/kinetics/KineticsFactory.h" #include "cantera/kinetics/GasKinetics.h" +#include "cantera/kinetics/Arrhenius.h" +#include "cantera/kinetics/ChebyshevRate.h" +#include "cantera/kinetics/Custom.h" +#include "cantera/kinetics/Falloff.h" +#include "cantera/kinetics/InterfaceRate.h" +#include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" #include "cantera/thermo/SurfPhase.h" -#include "cantera/kinetics/KineticsFactory.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/base/Array.h" @@ -35,9 +41,10 @@ TEST(Reaction, ElementaryFromYaml) " negative-A: true}"); auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_FALSE(R->usesThirdBody()); EXPECT_EQ(R->reactants.at("NO"), 1); EXPECT_EQ(R->products.at("N2"), 1); - EXPECT_EQ(R->type(), "reaction"); + EXPECT_EQ(R->type(), "Arrhenius"); EXPECT_EQ(R->rate()->type(), "Arrhenius"); EXPECT_FALSE(R->allow_negative_orders); @@ -47,7 +54,7 @@ TEST(Reaction, ElementaryFromYaml) EXPECT_DOUBLE_EQ(rate->activationEnergy(), 355 * 4184.0); } -TEST(Reaction, ThreeBodyFromYaml) +TEST(Reaction, ThreeBodyFromYaml1) { auto sol = newSolution("gri30.yaml", "", "None"); AnyMap rxn = AnyMap::fromYamlString( @@ -57,8 +64,28 @@ TEST(Reaction, ThreeBodyFromYaml) " efficiencies: {AR: 0.83, H2O: 5}}"); auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_TRUE(R->usesThirdBody()); + EXPECT_EQ(R->reactants.count("M"), (size_t) 0); + EXPECT_EQ(R->type(), "three-body-Arrhenius"); + EXPECT_DOUBLE_EQ(R->thirdBody()->efficiencies["H2O"], 5.0); + EXPECT_DOUBLE_EQ(R->thirdBody()->default_efficiency, 1.0); + + const auto& rate = std::dynamic_pointer_cast(R->rate()); + EXPECT_DOUBLE_EQ(rate->preExponentialFactor(), 1.2e11); +} + +TEST(Reaction, ThreeBodyFromYaml2) +{ + auto sol = newSolution("gri30.yaml", "", "None"); + AnyMap rxn = AnyMap::fromYamlString( + "{equation: 2 O + M = O2 + M," + " rate-constant: [1.20000E+17 cm^6/mol^2/s, -1, 0]," + " efficiencies: {AR: 0.83, H2O: 5}}"); + + auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_TRUE(R->usesThirdBody()); EXPECT_EQ(R->reactants.count("M"), (size_t) 0); - EXPECT_EQ(R->type(), "three-body"); + EXPECT_EQ(R->type(), "three-body-Arrhenius"); EXPECT_DOUBLE_EQ(R->thirdBody()->efficiencies["H2O"], 5.0); EXPECT_DOUBLE_EQ(R->thirdBody()->default_efficiency, 1.0); @@ -89,9 +116,10 @@ TEST(Reaction, FalloffFromYaml1) " efficiencies: {AR: 0.625}}"); auto R = newReaction(rxn, *(sol->kinetics())); - auto& FR = dynamic_cast(*R); - EXPECT_DOUBLE_EQ(FR.thirdBody()->efficiency("AR"), 0.625); - EXPECT_DOUBLE_EQ(FR.thirdBody()->efficiency("N2"), 1.0); + EXPECT_TRUE(R->usesThirdBody()); + EXPECT_EQ(R->type(), "falloff-SRI"); + EXPECT_DOUBLE_EQ(R->thirdBody()->efficiency("AR"), 0.625); + EXPECT_DOUBLE_EQ(R->thirdBody()->efficiency("N2"), 1.0); const auto rate = std::dynamic_pointer_cast(R->rate()); EXPECT_DOUBLE_EQ(rate->highRate().preExponentialFactor(), 7.91E+10); EXPECT_DOUBLE_EQ(rate->lowRate().preExponentialFactor(), 6.37E+14); @@ -110,9 +138,10 @@ TEST(Reaction, FalloffFromYaml2) " source: somewhere}"); auto R = newReaction(rxn, *(sol->kinetics())); - auto& FR = dynamic_cast(*R); - EXPECT_DOUBLE_EQ(FR.thirdBody()->efficiency("N2"), 1.0); - EXPECT_DOUBLE_EQ(FR.thirdBody()->efficiency("H2O"), 0.0); + EXPECT_TRUE(R->usesThirdBody()); + EXPECT_EQ(R->type(), "falloff-Troe"); + EXPECT_DOUBLE_EQ(R->thirdBody()->efficiency("N2"), 1.0); + EXPECT_DOUBLE_EQ(R->thirdBody()->efficiency("H2O"), 0.0); const auto rate = std::dynamic_pointer_cast(R->rate()); EXPECT_DOUBLE_EQ(rate->highRate().preExponentialFactor(), 6e11); EXPECT_DOUBLE_EQ(rate->lowRate().preExponentialFactor(), 1.04e20); @@ -138,9 +167,10 @@ TEST(Reaction, FalloffFromYaml3) " source: ARL-TR-5088}"); auto R = newReaction(rxn, *(sol->kinetics())); - auto& FR = dynamic_cast(*R); - EXPECT_DOUBLE_EQ(FR.thirdBody()->efficiency("N2"), 1.0); - EXPECT_DOUBLE_EQ(FR.thirdBody()->efficiency("H2O"), 5.0); + EXPECT_TRUE(R->usesThirdBody()); + EXPECT_EQ(R->type(), "falloff-Tsang"); + EXPECT_DOUBLE_EQ(R->thirdBody()->efficiency("N2"), 1.0); + EXPECT_DOUBLE_EQ(R->thirdBody()->efficiency("H2O"), 5.0); const auto rate = std::dynamic_pointer_cast(R->rate()); EXPECT_DOUBLE_EQ(rate->highRate().preExponentialFactor(), 8.3e17); EXPECT_DOUBLE_EQ(rate->lowRate().preExponentialFactor(), 3.57e26); @@ -164,6 +194,7 @@ TEST(Reaction, ChemicallyActivatedFromYaml) " low-P-rate-constant: [282320.078, 1.46878, -3270.56495]}"); auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_EQ(R->type(), "chemically-activated-Lindemann"); const auto& rate = std::dynamic_pointer_cast(R->rate()); EXPECT_DOUBLE_EQ(rate->highRate().preExponentialFactor(), 5.88e-14); EXPECT_DOUBLE_EQ(rate->lowRate().preExponentialFactor(), 2.82320078e2); @@ -184,6 +215,7 @@ TEST(Reaction, PlogFromYaml) "- {P: 1.01325 MPa, A: 1.680000e+16, b: -0.6, Ea: 14754.0}"); auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_FALSE(R->usesThirdBody()); const auto& rateMap = std::dynamic_pointer_cast(R->rate())->getRates(); std::vector> rates(rateMap.begin(), rateMap.end()); EXPECT_EQ(rates.size(), (size_t) 4); @@ -211,6 +243,7 @@ TEST(Reaction, ChebyshevFromYaml) " [-1.43220e-01, 7.71110e-02, 1.27080e-02, -6.41540e-04]]\n"); auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_FALSE(R->usesThirdBody()); EXPECT_EQ(R->reactants.size(), (size_t) 1); const auto& rate = std::dynamic_pointer_cast(R->rate()); double T = 1800; @@ -232,6 +265,7 @@ TEST(Reaction, BlowersMaselFromYaml) " negative-A: true}"); auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_FALSE(R->usesThirdBody()); EXPECT_EQ(R->reactants.at("H2"), 1); EXPECT_EQ(R->products.at("OH"), 1); @@ -265,10 +299,11 @@ TEST(Reaction, TwoTempPlasmaFromYaml) " rate-constant: [1.523e+27 cm^6/mol^2/s, -1.0, -100 K, 700 K]}"); auto R = newReaction(rxn, *(sol->kinetics())); - EXPECT_EQ(R->reactants.at("O2"), 2); + EXPECT_TRUE(R->usesThirdBody()); + EXPECT_EQ(R->thirdBody()->name(), "O2"); + EXPECT_EQ(R->reactants.at("O2"), 1); EXPECT_EQ(R->reactants.at("E"), 1); EXPECT_EQ(R->products.at("O2^-"), 1); - EXPECT_EQ(R->products.at("O2"), 1); const auto rate = std::dynamic_pointer_cast(R->rate()); EXPECT_DOUBLE_EQ(rate->preExponentialFactor(), 1.523e21); @@ -505,7 +540,7 @@ TEST_F(ReactionToYaml, elementary) soln = newSolution("h2o2.yaml", "", "None"); soln->thermo()->setState_TPY(1000, 2e5, "H2:1.0, O2:0.5, O:1e-8, OH:3e-8"); duplicateReaction(2); - EXPECT_EQ(duplicate->type(), "reaction"); + EXPECT_EQ(duplicate->type(), "Arrhenius"); EXPECT_EQ(duplicate->rate()->type(), "Arrhenius"); compareReactions(); } @@ -515,7 +550,7 @@ TEST_F(ReactionToYaml, threeBody) soln = newSolution("h2o2.yaml", "", "None"); soln->thermo()->setState_TPY(1000, 2e5, "H2:1.0, O2:0.5, O:1e-8, OH:3e-8, H:2e-7"); duplicateReaction(1); - EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + EXPECT_TRUE(std::dynamic_pointer_cast(duplicate->rate())); compareReactions(); } @@ -524,7 +559,9 @@ TEST_F(ReactionToYaml, TroeFalloff) soln = newSolution("h2o2.yaml", "", "None"); soln->thermo()->setState_TPY(1000, 2e5, "H2:1.0, O2:0.5, H2O2:1e-8, OH:3e-8"); duplicateReaction(21); - EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + auto rate = std::dynamic_pointer_cast(duplicate->rate()); + EXPECT_TRUE(rate); + EXPECT_FALSE(rate->chemicallyActivated()); compareReactions(); } @@ -533,7 +570,7 @@ TEST_F(ReactionToYaml, SriFalloff) soln = newSolution("sri-falloff.yaml"); soln->thermo()->setState_TPY(1000, 2e5, "R1A: 0.1, R1B:0.2, H: 0.2, R2:0.5"); duplicateReaction(0); - EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + EXPECT_TRUE(std::dynamic_pointer_cast(duplicate->rate())); compareReactions(); duplicateReaction(1); compareReactions(); @@ -544,7 +581,7 @@ TEST_F(ReactionToYaml, TsangFalloff) soln = newSolution("tsang-falloff.yaml"); soln->thermo()->setState_TPY(1000, 2e5, "NO:1.0, OH:1.0, H:1.0, CN:1.0"); duplicateReaction(0); - EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + EXPECT_TRUE(std::dynamic_pointer_cast(duplicate->rate())); compareReactions(); duplicateReaction(1); compareReactions(); @@ -555,7 +592,9 @@ TEST_F(ReactionToYaml, chemicallyActivated) soln = newSolution("chemically-activated-reaction.yaml"); soln->thermo()->setState_TPY(1000, 2e5, "H2:1.0, ch2o:0.1, ch3:1e-8, oh:3e-6"); duplicateReaction(0); - EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + auto rate = std::dynamic_pointer_cast(duplicate->rate()); + EXPECT_TRUE(rate); + EXPECT_TRUE(rate->chemicallyActivated()); compareReactions(); } @@ -635,12 +674,11 @@ TEST_F(ReactionToYaml, unconvertible2) TEST_F(ReactionToYaml, unconvertible3) { - FalloffReaction R( + Reaction R( {{"H2", 1}, {"OH", 1}}, {{"H2O", 1}, {"H", 1}}, - TroeRate( + shared_ptr(new TroeRate( ArrheniusRate(1e5, -1.0, 12.5), ArrheniusRate(1e5, -1.0, 12.5), - {0.562, 91.0, 5836.0, 8552.0}), - ThirdBody()); + {0.562, 91.0, 5836.0, 8552.0}))); AnyMap params = R.parameters(); UnitSystem U{"g", "cm", "mol"}; params.setUnits(U);