diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h index 471f3bb5109..aae4a96f313 100644 --- a/include/cantera/kinetics/BulkKinetics.h +++ b/include/cantera/kinetics/BulkKinetics.h @@ -11,12 +11,13 @@ #include "Kinetics.h" #include "RateCoeffMgr.h" +#include "ThirdBodyCalc.h" #include "cantera/kinetics/MultiRate.h" namespace Cantera { -class ElementaryReaction; +class ElementaryReaction2; //! Partial specialization of Kinetics for chemistry in a single bulk phase class BulkKinetics : public Kinetics @@ -45,9 +46,11 @@ class BulkKinetics : public Kinetics virtual void setMultiplier(size_t i, double f); virtual void invalidateCache(); + void addThirdBody(shared_ptr r); + protected: - virtual void addElementaryReaction(ElementaryReaction& r); - virtual void modifyElementaryReaction(size_t i, ElementaryReaction& rNew); + virtual void addElementaryReaction(ElementaryReaction2& r); + virtual void modifyElementaryReaction(size_t i, ElementaryReaction2& rNew); //! Vector of rate handlers std::vector> m_bulk_rates; @@ -62,12 +65,19 @@ class BulkKinetics : public Kinetics //! valued stoichiometries. vector_fp m_dn; - //! Activity concentrations, as calculated by - //! ThermoPhase::getActivityConcentrations + ThirdBodyCalc m_multi_concm; //!< used with MultiRate evaluator + vector_fp concm_multi_values; //!< concentrations of third-body collision partners + std::vector m_multi_indices; //!< reaction indices + + //! Third body concentrations + vector_fp m_concm; + + //! Activity concentrations, as calculated by ThermoPhase::getActivityConcentrations vector_fp m_act_conc; //! Physical concentrations, as calculated by ThermoPhase::getConcentrations vector_fp m_phys_conc; + vector_fp m_grt; bool m_ROP_ok; diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index 7f351a161a4..2ba5b768626 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -10,7 +10,6 @@ #define CT_GASKINETICS_H #include "BulkKinetics.h" -#include "ThirdBodyCalc.h" #include "FalloffMgr.h" #include "Reaction.h" @@ -87,7 +86,7 @@ class GasKinetics : public BulkKinetics ThirdBodyCalc m_falloff_concm; Rate1 m_plog_rates; - Rate1 m_cheb_rates; + Rate1 m_cheb_rates; Rate1 m_blowersmasel_rates; //! @name Reaction rate data @@ -107,16 +106,25 @@ class GasKinetics : public BulkKinetics void processFalloffReactions(); - void addThreeBodyReaction(ThreeBodyReaction& r); + // functions marked as deprecated below are only used for XML import and + // transitional reaction types are marked as '-legacy' + + //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach) + void addThreeBodyReaction(ThreeBodyReaction2& r); void addFalloffReaction(FalloffReaction& r); - void addPlogReaction(PlogReaction& r); - void addChebyshevReaction(ChebyshevReaction& r); + //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach) + void addPlogReaction(PlogReaction2& r); + //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach) + void addChebyshevReaction(ChebyshevReaction2& r); void addBlowersMaselReaction(BlowersMaselReaction& r); - void modifyThreeBodyReaction(size_t i, ThreeBodyReaction& r); + //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach) + void modifyThreeBodyReaction(size_t i, ThreeBodyReaction2& r); void modifyFalloffReaction(size_t i, FalloffReaction& r); - void modifyPlogReaction(size_t i, PlogReaction& r); - void modifyChebyshevReaction(size_t i, ChebyshevReaction& r); + //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach) + void modifyPlogReaction(size_t i, PlogReaction2& r); + //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach) + void modifyChebyshevReaction(size_t i, ChebyshevReaction2& r); void modifyBlowersMaselReaction(size_t i, BlowersMaselReaction& r); //! Update the equilibrium constants in molar units. diff --git a/include/cantera/kinetics/MultiRate.h b/include/cantera/kinetics/MultiRate.h index e2fad2002d1..0bdb00c5962 100644 --- a/include/cantera/kinetics/MultiRate.h +++ b/include/cantera/kinetics/MultiRate.h @@ -1,8 +1,5 @@ /** * @file MultiRate.h - * - * @warning This file is an experimental part of the %Cantera API and - * may be changed or removed without notice. */ // This file is part of Cantera. See License.txt in the top-level directory or @@ -26,9 +23,6 @@ namespace Cantera * which can be updated using information of a single `ThermoPhase`. * `InterfaceKinetics` will require access to an entire `Kinetics` object * or the underlying `vector` vector (e.g. `m_thermo`). - * - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. */ class MultiRateBase { @@ -47,37 +41,40 @@ class MultiRateBase virtual bool replace(const size_t rxn_index, ReactionRateBase& rate) = 0; + //! Update number of species + //! @param n_species number of species + virtual void resizeSpecies(size_t n_species) = 0; + //! Evaluate all rate constants handled by the evaluator //! @param bulk object representing bulk phase //! @param kf array of rate constants virtual void getRateConstants(const ThermoPhase& bulk, - double* kf) const = 0; + double* kf, double* concm) const = 0; //! Update data common to reaction rates of a specific type //! @param bulk object representing bulk phase - virtual void update(const ThermoPhase& bulk) = 0; + //! @param concm effective third-body concentrations + //! @TODO enable more generic handling of non-trivial concentration dependencies + virtual void update(const ThermoPhase& bulk, double* concm) = 0; }; //! A class template handling all reaction rates specific to `BulkKinetics`. -/** - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. - */ template class MultiBulkRates final : public MultiRateBase { public: virtual void add(const size_t rxn_index, - ReactionRateBase& rate) override { - m_indices[rxn_index] = m_rates.size(); - m_rates.push_back(dynamic_cast(rate)); - m_rxn.push_back(rxn_index); + ReactionRateBase& rate) override + { + m_indices[rxn_index] = m_rxn_rates.size(); + m_rxn_rates.emplace_back(rxn_index, dynamic_cast(rate)); } virtual bool replace(const size_t rxn_index, - ReactionRateBase& rate) override { - if (!m_rates.size()) { + ReactionRateBase& rate) override + { + if (!m_rxn_rates.size()) { throw CanteraError("MultiBulkRate::replace", "Invalid operation: cannot replace rate object " "in empty rate handler."); @@ -85,39 +82,46 @@ class MultiBulkRates final : public MultiRateBase throw CanteraError("MultiBulkRate::replace", "Invalid operation: cannot replace rate object of type '{}' " "with a new rate of type '{}'.", - m_rates[0].type(), rate.type()); + m_rxn_rates.at(0).second.type(), rate.type()); } if (m_indices.find(rxn_index) != m_indices.end()) { size_t j = m_indices[rxn_index]; - m_rates[j] = dynamic_cast(rate); + m_rxn_rates.at(j).second = dynamic_cast(rate); return true; } return false; } + virtual void resizeSpecies(size_t n_species) + { + m_shared.resizeSpecies(n_species); + } + virtual void getRateConstants(const ThermoPhase& bulk, - double* kf) const override { - for (size_t i = 0; i < m_rates.size(); i++) { - kf[m_rxn[i]] = m_rates[i].eval(m_shared); + double* kf, double* concm) const override + { + for (const auto& rxn : m_rxn_rates) { + kf[rxn.first] = rxn.second.eval(m_shared, concm[rxn.first]); } } - virtual void update(const ThermoPhase& bulk) override { + virtual void update(const ThermoPhase& bulk, double* concm) override + { // update common data once for each reaction type m_shared.update(bulk); - if (RateType::uses_update()) { + if (RateType::usesUpdate()) { // update reaction-specific data for each reaction. This loop // is efficient as all function calls are de-virtualized, and // all of the rate objects are contiguous in memory - for (auto& rate : m_rates) { - rate.update(m_shared, bulk); + for (auto& rxn : m_rxn_rates) { + rxn.second.update(m_shared, concm[rxn.first]); } } } protected: - std::vector m_rates; //! Reaction rate objects - std::vector m_rxn; //! Index within overall rate vector + //! Vector of pairs of reaction rates indices and reaction rates + std::vector> m_rxn_rates; std::map m_indices; //! Mapping of indices DataType m_shared; }; diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 6922d4930b9..84f6511c0b1 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -19,6 +19,7 @@ namespace Cantera class Kinetics; class Falloff; class XML_Node; +class ThirdBody; //! Abstract base class which stores data about a reaction and its rate //! parameterization so that it can be added to a Kinetics object. @@ -65,6 +66,9 @@ class Reaction //! contained in the #input attribute. AnyMap parameters(bool withInput=true) const; + //! Set up reaction based on AnyMap *node* + virtual void setParameters(const AnyMap& node, const Kinetics& kin); + //! Get validity flag of reaction bool valid() const { return m_valid; @@ -132,6 +136,26 @@ class Reaction //! where the reaction occurs. Units rate_units; + //! Get reaction rate pointer + shared_ptr rate() { + return m_rate; + } + + //! Set reaction rate pointer + void setRate(shared_ptr rate) { + m_rate = rate; + } + + //! Get pointer to third-body + shared_ptr thirdBody() { + return m_third_body; + } + + //! Indicate whether object uses legacy framework + bool usesLegacy() const { + return !m_rate; + } + protected: //! Store the parameters of a Reaction needed to reconstruct an identical //! object using the newReaction(AnyMap&, Kinetics&) function. Does not @@ -148,59 +172,47 @@ class Reaction //! @param kin Kinetics object virtual std::pair, bool> undeclaredThirdBodies(const Kinetics& kin) const; -}; - -//! An intermediate class used to avoid naming conflicts of 'rate' member -//! variables and getters (see `ElementaryReaction`, `PlogReaction` and -//! `ChebyshevReaction`). -class Reaction2 : public Reaction -{ -public: - //! Get reaction rate pointer - shared_ptr rate() { - return m_rate; - } - - //! Set reaction rate pointer - void setRate(shared_ptr rate) { - m_rate = rate; - } - - //! Set up reaction based on AnyMap node - virtual void setParameters(const AnyMap& node, const Kinetics& kin); - -protected: //! Reaction rate used by generic reactions shared_ptr m_rate; + + //! Relative efficiencies of third-body species in enhancing the reaction + //! rate (if applicable) + shared_ptr m_third_body; }; //! A reaction which follows mass-action kinetics with a modified Arrhenius //! reaction rate. -class ElementaryReaction : public Reaction +class ElementaryReaction2 : public Reaction { public: - ElementaryReaction(); - ElementaryReaction(const Composition& reactants, const Composition products, - const Arrhenius& rate); + ElementaryReaction2(); + ElementaryReaction2(const Composition& reactants, const Composition products, + const Arrhenius& rate); virtual void validate(); virtual void getParameters(AnyMap& reactionNode) const; virtual std::string type() const { - return "elementary"; + return "elementary-legacy"; } Arrhenius rate; bool allow_negative_pre_exponential_factor; }; + //! A class for managing third-body efficiencies, including default values class ThirdBody { public: explicit ThirdBody(double default_efficiency=1.0); + ThirdBody(const AnyMap& node); + + //! Set third-body efficiencies from AnyMap *node* + void setEfficiencies(const AnyMap& node); + //! Get the third-body efficiency for species *k* double efficiency(const std::string& k) const; @@ -212,17 +224,18 @@ class ThirdBody double default_efficiency; }; + //! A reaction with a non-reacting third body "M" that acts to add or remove //! energy from the reacting species -class ThreeBodyReaction : public ElementaryReaction +class ThreeBodyReaction2 : public ElementaryReaction2 { public: - ThreeBodyReaction(); - ThreeBodyReaction(const Composition& reactants, const Composition& products, - const Arrhenius& rate, const ThirdBody& tbody); + ThreeBodyReaction2(); + ThreeBodyReaction2(const Composition& reactants, const Composition& products, + const Arrhenius& rate, const ThirdBody& tbody); virtual std::string type() const { - return "three-body"; + return "three-body-legacy"; } virtual std::string reactantString() const; @@ -241,6 +254,7 @@ class ThreeBodyReaction : public ElementaryReaction undeclaredThirdBodies(const Kinetics& kin) const; }; + //! A reaction that is first-order in [M] at low pressure, like a third-body //! reaction, but zeroth-order in [M] as pressure increases. class FalloffReaction : public Reaction @@ -285,6 +299,7 @@ class FalloffReaction : public Reaction const Kinetics& kin) const; }; + //! A reaction where the rate decreases as pressure increases due to collisional //! stabilization of a reaction intermediate. Like a FalloffReaction, except //! that the forward rate constant is written as being proportional to the low- @@ -305,17 +320,18 @@ class ChemicallyActivatedReaction : public FalloffReaction virtual void getParameters(AnyMap& reactionNode) const; }; + //! A pressure-dependent reaction parameterized by logarithmically interpolating //! between Arrhenius rate expressions at various pressures. -class PlogReaction : public Reaction +class PlogReaction2 : public Reaction { public: - PlogReaction(); - PlogReaction(const Composition& reactants, const Composition& products, - const Plog& rate); + PlogReaction2(); + PlogReaction2(const Composition& reactants, const Composition& products, + const Plog& rate); virtual std::string type() const { - return "pressure-dependent-Arrhenius"; + return "pressure-dependent-Arrhenius-legacy"; } virtual void validate(); @@ -324,62 +340,22 @@ class PlogReaction : public Reaction Plog rate; }; + //! A pressure-dependent reaction parameterized by a bi-variate Chebyshev //! polynomial in temperature and pressure -class ChebyshevReaction : public Reaction +class ChebyshevReaction2 : public Reaction { public: - ChebyshevReaction(); - ChebyshevReaction(const Composition& reactants, const Composition& products, - const ChebyshevRate& rate); + ChebyshevReaction2(); + ChebyshevReaction2(const Composition& reactants, const Composition& products, + const Chebyshev& rate); virtual void getParameters(AnyMap& reactionNode) const; virtual std::string type() const { - return "Chebyshev"; + return "Chebyshev-legacy"; } - ChebyshevRate rate; -}; - -//! A reaction which follows mass-action kinetics with a custom reaction rate -//! defined in Python. -/** - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. - */ -class CustomFunc1Reaction : public Reaction2 -{ -public: - CustomFunc1Reaction(); - - virtual void setParameters(const AnyMap& node, const Kinetics& kin); - - virtual std::string type() const { - return "custom-rate-function"; - } -}; - - -//! A reaction which follows mass-action kinetics with a modified Arrhenius -//! reaction rate. -/** - * Alternative elementary reaction based on ReactionRate. - * - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. - */ -class TestReaction : public Reaction2 -{ -public: - TestReaction(); - - virtual std::string type() const { - return "elementary-new"; - } - - virtual void setParameters(const AnyMap& node, const Kinetics& kin); - - bool allow_negative_pre_exponential_factor; + Chebyshev rate; }; @@ -398,8 +374,9 @@ struct CoverageDependency double m; //!< exponent for power law dependence on coverage [dimensionless] }; + //! A reaction occurring on an interface (i.e. a SurfPhase or an EdgePhase) -class InterfaceReaction : public ElementaryReaction +class InterfaceReaction : public ElementaryReaction2 { public: InterfaceReaction(); @@ -433,6 +410,7 @@ class InterfaceReaction : public ElementaryReaction std::string sticking_species; }; + //! An interface reaction which involves charged species class ElectrochemicalReaction : public InterfaceReaction { @@ -448,6 +426,7 @@ class ElectrochemicalReaction : public InterfaceReaction bool exchange_current_density_formulation; }; + //! A reaction with rate parameters for Blowers-Masel approximation class BlowersMaselReaction: public Reaction { @@ -467,6 +446,7 @@ class BlowersMaselReaction: public Reaction bool allow_negative_pre_exponential_factor; }; + //! A reaction occurring on an interface (i.e. a SurfPhase or an EdgePhase) //! with the rate calculated with Blowers-Masel approximation. class BlowersMaselInterfaceReaction : public BlowersMaselReaction @@ -502,6 +482,123 @@ class BlowersMaselInterfaceReaction : public BlowersMaselReaction std::string sticking_species; }; + +//! A reaction which follows mass-action kinetics with a modified Arrhenius +//! reaction rate. +class ElementaryReaction3 : public Reaction +{ +public: + ElementaryReaction3(); + ElementaryReaction3(const Composition& reactants, const Composition& products, + const ArrheniusRate& rate); + + ElementaryReaction3(const AnyMap& node, const Kinetics& kin); + + virtual std::string type() const { + return "elementary"; + } + virtual void getParameters(AnyMap& reactionNode) const; +}; + + +//! A reaction with a non-reacting third body "M" that acts to add or remove +//! energy from the reacting species +class ThreeBodyReaction3 : public ElementaryReaction3 +{ +public: + ThreeBodyReaction3(); + ThreeBodyReaction3(const Composition& reactants, const Composition& products, + const ArrheniusRate& rate, const ThirdBody& tbody); + + ThreeBodyReaction3(const AnyMap& node, const Kinetics& kin); + + virtual std::string type() const { + return "three-body"; + } + + bool detectEfficiencies(); + virtual void calculateRateCoeffUnits(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; + + bool specified_collision_partner = false; //!< Input specifies collision partner +}; + + +//! A pressure-dependent reaction parameterized by logarithmically interpolating +//! between Arrhenius rate expressions at various pressures. +class PlogReaction3 : public Reaction +{ +public: + PlogReaction3(); + PlogReaction3(const Composition& reactants, const Composition& products, + const PlogRate& rate); + + PlogReaction3(const AnyMap& node, const Kinetics& kin); + + virtual std::string type() const { + return "pressure-dependent-Arrhenius"; + } + virtual void getParameters(AnyMap& reactionNode) const; +}; + + +//! A pressure-dependent reaction parameterized by a bi-variate Chebyshev +//! polynomial in temperature and pressure +class ChebyshevReaction3 : public Reaction +{ +public: + ChebyshevReaction3(); + ChebyshevReaction3(const Composition& reactants, const Composition& products, + const ChebyshevRate3& rate); + + ChebyshevReaction3(const AnyMap& node, const Kinetics& kin); + + virtual std::string type() const { + return "Chebyshev"; + } + + virtual void getParameters(AnyMap& reactionNode) const; +}; + + +//! A reaction which follows mass-action kinetics with a custom reaction rate +//! defined in Python. +/** + * @warning This class is an experimental part of the %Cantera API and + * may be changed or removed without notice. + */ +class CustomFunc1Reaction : public Reaction +{ +public: + CustomFunc1Reaction(); + CustomFunc1Reaction(const Composition& reactants, const Composition& products, + const CustomFunc1Rate& rate); + + CustomFunc1Reaction(const AnyMap& node, const Kinetics& kin); + + virtual std::string type() const { + return "custom-rate-function"; + } +}; + + +#ifdef CT_NO_LEGACY_REACTIONS_26 +typedef ElementaryReaction3 ElementaryReaction; +typedef ThreeBodyReaction3 ThreeBodyReaction; +typedef PlogReaction3 PlogReaction; +typedef ChebyshevReaction3 ChebyshevReaction; +#else +typedef ElementaryReaction2 ElementaryReaction; +typedef ThreeBodyReaction2 ThreeBodyReaction; +typedef PlogReaction2 PlogReaction; +typedef ChebyshevReaction2 ChebyshevReaction; +#endif + + //! Create Reaction objects for all `` nodes in an XML document. //! //! The `` nodes are assumed to be children of the `` @@ -538,14 +635,14 @@ void parseReactionEquation(Reaction& R, const AnyValue& equation, // declarations of setup functions void setupReaction(Reaction& R, const XML_Node& rxn_node); -void setupElementaryReaction(ElementaryReaction&, const XML_Node&); +void setupElementaryReaction(ElementaryReaction2&, const XML_Node&); //! @internal May be changed without notice in future versions -void setupElementaryReaction(ElementaryReaction&, const AnyMap&, +void setupElementaryReaction(ElementaryReaction2&, const AnyMap&, const Kinetics&); -void setupThreeBodyReaction(ThreeBodyReaction&, const XML_Node&); -//! @internal May be changed without notice in future versions -void setupThreeBodyReaction(ThreeBodyReaction&, const AnyMap&, +void setupThreeBodyReaction(ThreeBodyReaction2&, const XML_Node&); +//! @deprecated Cantera 2.6 (replaced by setParameters) +void setupThreeBodyReaction(ThreeBodyReaction2&, const AnyMap&, const Kinetics&); void setupFalloffReaction(FalloffReaction&, const XML_Node&); @@ -556,13 +653,13 @@ void setupFalloffReaction(FalloffReaction&, const AnyMap&, void setupChemicallyActivatedReaction(ChemicallyActivatedReaction&, const XML_Node&); -void setupPlogReaction(PlogReaction&, const XML_Node&); -//! @internal May be changed without notice in future versions -void setupPlogReaction(PlogReaction&, const AnyMap&, const Kinetics&); +void setupPlogReaction(PlogReaction2&, const XML_Node&); +//! @deprecated Cantera 2.6 (replaced by setParameters) +void setupPlogReaction(PlogReaction2&, const AnyMap&, const Kinetics&); -void setupChebyshevReaction(ChebyshevReaction&, const XML_Node&); -//! @internal May be changed without notice in future versions -void setupChebyshevReaction(ChebyshevReaction&, const AnyMap&, +void setupChebyshevReaction(ChebyshevReaction2&, const XML_Node&); +//! @deprecated Cantera 2.6 (replaced by setParameters) +void setupChebyshevReaction(ChebyshevReaction2&, const AnyMap&, const Kinetics&); void setupInterfaceReaction(InterfaceReaction&, const XML_Node&); diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h index f5ae68317f4..b7117c8ddd3 100644 --- a/include/cantera/kinetics/ReactionData.h +++ b/include/cantera/kinetics/ReactionData.h @@ -1,8 +1,5 @@ /** * @file ReactionData.h - * - * @warning This file is an experimental part of the %Cantera API and - * may be changed or removed without notice. */ // This file is part of Cantera. See License.txt in the top-level directory or @@ -23,59 +20,117 @@ class ThermoPhase; /** * The data container `ArrheniusData` holds precalculated data common to * all `ArrheniusRate` objects. - * - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. */ struct ArrheniusData { ArrheniusData() : m_temperature(1.), m_logT(0.), m_recipT(1.) {} - //! Constructor based on temperature *T* - ArrheniusData(double T) { update(T); } + //! Update data container based on temperature *T* + void update(double T) + { + m_temperature = T; + m_logT = std::log(T); + m_recipT = 1./T; + } + + //! Update data container based on temperature *T* and pressure *P* + void update(double T, double P) { update(T); } + + //! Update data container based on *bulk* phase state + void update(const ThermoPhase& bulk); + + //! Update number of species; unused + void resizeSpecies(size_t n_species) {} + + double m_temperature; //!< temperature + double m_logT; //!< logarithm of temperature + double m_recipT; //!< inverse of temperature +}; - //! Constructor based on temperature *T* and pressure *P* - ArrheniusData(double T, double P) { update(T); }; - //! Constructor accessing *bulk* phase definitions - ArrheniusData(const ThermoPhase& bulk) { update(bulk); } +//! Data container holding shared data specific to PlogRate +/** + * The data container `PlogData` holds precalculated data common to + * all `PlogRate` objects. + */ +struct PlogData +{ + PlogData() : m_temperature(1.), m_logT(0.), m_recipT(1.), m_logP(0.) {} + + //! Update data container based on temperature *T* (raises exception) + void update(double T); - void update(double T) { + //! Update data container based on temperature *T* and *P* + void update(double T, double P) + { m_temperature = T; m_logT = std::log(T); m_recipT = 1./T; - } + m_logP = std::log(P); + } + //! Update data container based on *bulk* phase state void update(const ThermoPhase& bulk); + //! Update number of species; unused + void resizeSpecies(size_t n_species) {} + double m_temperature; //!< temperature double m_logT; //!< logarithm of temperature - double m_recipT; //!< inverse of temperature + double m_recipT; //!< inverse of temperature + double m_logP; //!< logarithm of pressure }; -//! Data container holding shared data specific to CustomFunc1Rate +//! Data container holding shared data specific to ChebyshevRate /** - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. + * The data container `ChebyshevData` holds precalculated data common to + * all `ChebyshevRate3` objects. */ -struct CustomFunc1Data +struct ChebyshevData { - CustomFunc1Data() : m_temperature(1.) {} + ChebyshevData() : m_temperature(1.), m_recipT(1.), m_log10P(0.) {} - //! Constructor based on temperature *T* - CustomFunc1Data(double T) { update(T); } + //! Update data container based on temperature *T* (raises exception) + void update(double T); + + //! Update data container based on temperature *T* and *P* + void update(double T, double P) + { + m_temperature = T; + m_recipT = 1./T; + m_log10P = std::log10(P); + } - //! Constructor based on temperature *T* and pressure *P* - CustomFunc1Data(double T, double P) { update(T); }; + //! Update data container based on *bulk* phase state + void update(const ThermoPhase& bulk); - //! Constructor accessing *bulk* phase definitions - CustomFunc1Data(const ThermoPhase& bulk) { update(bulk); } + //! Update number of species; unused + void resizeSpecies(size_t n_species) {} + double m_temperature; //!< temperature + double m_recipT; //!< inverse of temperature + double m_log10P; //!< base 10 logarithm of pressure +}; + + +//! Data container holding shared data specific to CustomFunc1Rate +struct CustomFunc1Data +{ + CustomFunc1Data() : m_temperature(1.) {} + + //! Update data container based on temperature *T* void update(double T) { m_temperature = T; } + //! Update data container based on temperature *T* and pressure *P* + void update(double T, double P) { update(T); } + + //! Update data container based on *bulk* phase state void update(const ThermoPhase& bulk); + //! Update number of species; unused + void resizeSpecies(size_t n_species) {} + double m_temperature; //!< temperature }; diff --git a/include/cantera/kinetics/ReactionFactory.h b/include/cantera/kinetics/ReactionFactory.h index 6cf0ff7e09c..45a956fae6a 100644 --- a/include/cantera/kinetics/ReactionFactory.h +++ b/include/cantera/kinetics/ReactionFactory.h @@ -26,13 +26,8 @@ namespace Cantera * @endcode * * @ingroup reactionGroup - * - * @todo Once support for XML is phased out, a change of the base class type - * to `Factory` will ensure that - * `reg` and `create` methods provided by the base factory class can be - * used. Thus, separate `setup_AnyMap` and other methods can be avoided. */ -class ReactionFactory : public Factory +class ReactionFactory : public Factory { public: /** @@ -54,73 +49,51 @@ class ReactionFactory : public Factory s_factory = 0; } - //! Set up reaction based on XML node information - /** - * Kept separate from base `create` method as the factory has to handle both - * XML and `AnyMap` input until support for the former is removed. - */ - void setup_XML(std::string name, - Reaction* R, const XML_Node& rxn_node) { - try { - m_xml_setup.at(name)(R, rxn_node); - } catch (std::out_of_range&) { - throw CanteraError("ReactionFactory::setup_XML", - "No such type: '{}'", name); - } - } +private: + //! Pointer to the single instance of the factory + static ReactionFactory* s_factory; - //! Register a new object initializer function (from XML node) - /** - * Kept separate from base `reg` method as the factory has to handle both - * XML and `AnyMap` input until support for the former is removed. - */ - void reg_XML(const std::string& name, - std::function f) { - m_xml_setup[name] = f; - } + //! default constructor, which is defined as private + ReactionFactory(); - //! Set up reaction based on AnyMap node information + //! Mutex for use when calling the factory + static std::mutex reaction_mutex; +}; + +class ReactionFactoryXML : public Factory +{ +public: /** - * @todo call setup directly from constructor after removal of XML support. + * 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. */ - void setup_AnyMap(std::string name, - Reaction* R, const AnyMap& node, const Kinetics& kin) { - try { - m_anymap_setup.at(name)(R, node, kin); - } catch (std::out_of_range&) { - throw CanteraError("ReactionFactory::setup_AnyMap", - "No such type: '{}'", name); + static ReactionFactoryXML* factory() { + std::unique_lock lock(reaction_mutex); + if (!s_factory) { + s_factory = new ReactionFactoryXML; } + return s_factory; } - //! Register a new object initializer function (from AnyMap node) - /** - * @todo can be handled by `reg` after removal of XML support. - */ - void reg_AnyMap(const std::string& name, - std::function f) { - m_anymap_setup[name] = f; + 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; + static ReactionFactoryXML* s_factory; //! default constructor, which is defined as private - ReactionFactory(); + ReactionFactoryXML(); //! Mutex for use when calling the factory static std::mutex reaction_mutex; - - //! map of XML initializers - std::unordered_map> m_xml_setup; - - //! map of AnyMap initializers - std::unordered_map> m_anymap_setup; }; + //! Create a new empty Reaction object /*! * @param type string identifying type of reaction. diff --git a/include/cantera/kinetics/ReactionRate.h b/include/cantera/kinetics/ReactionRate.h index 52635465b84..53448d84cd8 100644 --- a/include/cantera/kinetics/ReactionRate.h +++ b/include/cantera/kinetics/ReactionRate.h @@ -11,6 +11,8 @@ #ifndef CT_REACTIONRATE_H #define CT_REACTIONRATE_H +#include "cantera/base/AnyMap.h" +#include "cantera/base/Units.h" #include "cantera/kinetics/RxnRates.h" #include "cantera/kinetics/ReactionData.h" #include "cantera/base/ctexceptions.h" @@ -19,8 +21,6 @@ namespace Cantera { class Func1; -class AnyMap; - //! Abstract base class for reaction rate definitions /** @@ -32,18 +32,32 @@ class AnyMap; * Methods defined for the abstract base class are not aware of specialized * data handlers defined by the template class `ReactionRate` * and thus can be exposed to the API. - * - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. */ class ReactionRateBase { public: + ReactionRateBase() : units(0.) {} virtual ~ReactionRateBase() {} +public: //! Identifier of reaction type virtual std::string type() const = 0; + //! Update reaction rate data based on temperature + //! @param T temperature [K] + virtual void update(double T) = 0; + + //! Update reaction rate data based on temperature and pressure + //! @param T temperature [K] + //! @param P pressure [Pa] + //! @param concm third-body concentration (if applicable) + virtual void update(double T, double P, double concm=0.) = 0; + + //! Update reaction rate data based on bulk phase + //! @param bulk object representing bulk phase + //! @param concm third-body concentration (if applicable) + virtual void update(const ThermoPhase& bulk, double concm=0.) = 0; + //! Evaluate reaction rate based on temperature //! @param T temperature [K] virtual double eval(double T) const = 0; @@ -51,11 +65,13 @@ class ReactionRateBase //! Evaluate reaction rate based on temperature and pressure //! @param T temperature [K] //! @param P pressure [Pa] - virtual double eval(double T, double P) const = 0; + //! @param concm third-body concentration (if applicable) + virtual double eval(double T, double P, double concm=0.) const = 0; //! Evaluate reaction rate based on bulk phase //! @param bulk object representing bulk phase - virtual double eval(const ThermoPhase& bulk) const = 0; + //! @param concm third-body concentration (if applicable) + virtual double eval(const ThermoPhase& bulk, double concm=0.) const = 0; //! Evaluate reaction rate derivative based on temperature //! @param T temperature [K] @@ -68,7 +84,42 @@ class ReactionRateBase //! Evaluate reaction rate derivative based on bulk phase //! @param bulk object representing bulk phase - virtual double ddT(const ThermoPhase& bulk) const = 0; + //! @param concm third-body concentration (if applicable) + virtual double ddT(const ThermoPhase& bulk, double concm=0.) const = 0; + + //! Validate the reaction rate expression + virtual void validate(const std::string& equation) = 0; + + //! Return the parameters such that an identical Reaction could be reconstructed + //! using the newReaction() function. Behavior specific to derived classes is + //! handled by the getParameters() method. + //! @param rate_units units used for rate parameters + AnyMap parameters(const Units& rate_units) const; + + //! Return parameters using original unit system + AnyMap parameters() const; + + //! Set parameters + //! @param node AnyMap object containing reaction rate specification + //! @param rate_units Description of units used for rate parameters + virtual void setParameters(const AnyMap& node, const Units& rate_units); + +protected: + //! Get parameters + //! Store the parameters of a ReactionRate needed to reconstruct an identical + //! object. Does not include user-defined fields available in the #input map. + virtual void getParameters(AnyMap& rateNode, const Units& rate_units) const { + throw CanteraError("ReactionRate::getParameters", + "Not implemented by derived ReactionRate object."); + } + + //! Input data used for specific models + AnyMap input; + + //! The units of the rate constant. These are determined by the units of the + //! standard concentration of the reactant species' phases of the phase + //! where the reaction occurs. + Units units; }; @@ -76,9 +127,6 @@ class ReactionRateBase /** * This class template ensures that derived objects are aware of specialized * data types, which are passed by MultiRateBase evaluators. - * - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. */ template class ReactionRate : public ReactionRateBase @@ -88,77 +136,282 @@ class ReactionRate : public ReactionRateBase //! Update information specific to reaction //! @param shared_data data shared by all reactions of a given type - virtual void update(const DataType& shared_data, const ThermoPhase& bulk) {} + virtual void update(const DataType& shared_data, + double concm=0.) {} + + virtual void update(double T) override { + DataType data; + data.update(T); + update(data); + } + + virtual void update(double T, double P, double concm=0.) override { + DataType data; + data.update(T, P); + update(data, concm); + } + + virtual void update(const ThermoPhase& bulk, double concm=0.) override { + DataType data; + data.update(bulk); + update(data); + } //! Evaluate reaction rate //! @param shared_data data shared by all reactions of a given type - virtual double eval(const DataType& shared_data) const = 0; + //! @param concm third-body concentration (if applicable) + virtual double eval(const DataType& shared_data, double concm=0.) const = 0; virtual double eval(double T) const override { - return eval(DataType(T)); + DataType data; + data.update(T); + return eval(data); } - virtual double eval(double T, double P) const override { - return eval(DataType(T, P)); + virtual double eval(double T, double P, double concm=0.) const override { + DataType data; + data.update(T, P); + return eval(data, concm); } - virtual double eval(const ThermoPhase& bulk) const override { - return eval(DataType(bulk)); + virtual double eval(const ThermoPhase& bulk, double concm=0.) const override { + DataType data; + data.update(bulk); + return eval(data, concm); } //! Evaluate derivative of reaction rate with respect to temperature //! @param shared_data data shared by all reactions of a given type - virtual double ddT(const DataType& shared_data) const { + virtual double ddT(const DataType& shared_data, double concm=0.) const { throw CanteraError("ReactionRate::ddT", "Not implemented by derived ReactionRate object."); } virtual double ddT(double T) const override { - return ddT(DataType(T)); + DataType data; + data.update(T); + return ddT(data); } virtual double ddT(double T, double P) const override { - return ddT(DataType(T, P)); + DataType data; + data.update(T, P); + return ddT(data); } - virtual double ddT(const ThermoPhase& bulk) const override { - return ddT(DataType(bulk)); + virtual double ddT(const ThermoPhase& bulk, double concm=0.) const override { + DataType data; + data.update(bulk); + return ddT(data, concm); } }; -//! Arrhenius reaction rate type depends only on temperature +//! The ArrheniusRate reaction rate type depends only on temperature /** - * Wrapped Arrhenius rate. + * A reaction rate coefficient of the following form. * - * @warning This class is an experimental part of the %Cantera API and - * may be changed or removed without notice. + * \f[ + * k_f = A T^b \exp (-E/RT) + * \f] + * + * Note: ArrheniusRate acts as a wrapper for the Arrhenius class. */ class ArrheniusRate final : public ReactionRate, public Arrhenius { public: + //! Default constructor. ArrheniusRate(); + //! Constructor. + //! @param A pre-exponential. The unit system is + //! (kmol, m, s). The actual units depend on the reaction + //! order and the dimensionality (surface or bulk). + //! @param b Temperature exponent. Non-dimensional. + //! @param E Activation energy. J/kmol. ArrheniusRate(double A, double b, double E); - //! Constructor - //! @param node AnyMap object containing reaction rate specification - //! @param rate_units Description of units used for rate parameters + //! Constructor using AnyMap content + //! @param node AnyMap containing rate information + //! @param rate_units unit definitions used for rate information ArrheniusRate(const AnyMap& node, const Units& rate_units); + //! Constructor using AnyMap content + //! @param node AnyMap containing rate information + ArrheniusRate(const AnyMap& node); + + //! Constructor based on Arrhenius object + //! @param arr Arrhenius object + //! @param allow_negative_A allow negative pre-exponential factor + //! (optional, default is false) + ArrheniusRate(const Arrhenius& arr, bool allow_negative_A=false); + + virtual void setParameters(const AnyMap& node, const Units& rate_units) override; + virtual void getParameters(AnyMap& rateNode, + const Units& rate_units) const override; + virtual std::string type() const override { return "ArrheniusRate"; } //! Update information specific to reaction - static bool uses_update() { return false; } + static bool usesUpdate() { return false; } - virtual double eval(const ArrheniusData& shared_data) const override { + virtual double eval(const ArrheniusData& shared_data, + double concm=0.) const override { return updateRC(shared_data.m_logT, shared_data.m_recipT); } - virtual double ddT(const ArrheniusData& shared_data) const override { + virtual double ddT(const ArrheniusData& shared_data, + double concm=0.) const override { return updateRC(shared_data.m_logT, shared_data.m_recipT) * (m_b + m_E * shared_data.m_recipT) * shared_data.m_recipT; } + + //! Return the activation energy [J/kmol] + double activationEnergy() const { + return m_E * GasConstant; + } + + virtual void validate(const std::string& equation) override; + + bool allow_negative_pre_exponential_factor; +}; + + +//! Pressure-dependent reaction rate expressed by logarithmically interpolating +//! between Arrhenius rate expressions at various pressures. +/*! + * Given two rate expressions at two specific pressures: + * + * * \f$ P_1: k_1(T) = A_1 T^{b_1} e^{-E_1 / RT} \f$ + * * \f$ P_2: k_2(T) = A_2 T^{b_2} e^{-E_2 / RT} \f$ + * + * The rate at an intermediate pressure \f$ P_1 < P < P_2 \f$ is computed as + * \f[ + * \log k(T,P) = \log k_1(T) + \bigl(\log k_2(T) - \log k_1(T)\bigr) + * \frac{\log P - \log P_1}{\log P_2 - \log P_1} + * \f] + * Multiple rate expressions may be given at the same pressure, in which case + * the rate used in the interpolation formula is the sum of all the rates given + * at that pressure. For pressures outside the given range, the rate expression + * at the nearest pressure is used. + */ +class PlogRate final : public ReactionRate, public Plog +{ +public: + PlogRate() {} + + //! Constructor from Arrhenius rate expressions at a set of pressures + explicit PlogRate(const std::multimap& rates); + + //! Constructor using AnyMap content + //! @param node AnyMap containing rate information + //! @param rate_units unit definitions used for rate information + PlogRate(const AnyMap& node, const Units& rate_units); + + //! Constructor using AnyMap content + //! @param node AnyMap containing rate information + PlogRate(const AnyMap& node); + + virtual std::string type() const override { return "PlogRate"; } + + virtual void setParameters(const AnyMap& node, const Units& rate_units) override; + virtual void getParameters(AnyMap& rateNode, + const Units& rate_units) const override; + + //! Update information specific to reaction + static bool usesUpdate() { return true; } + + virtual void update(const PlogData& shared_data, double concm=0.) override { + update_C(&shared_data.m_logP); + } + + virtual double eval(const PlogData& shared_data, + double concm=0.) const override { + return updateRC(shared_data.m_logT, shared_data.m_recipT); + } + + virtual void validate(const std::string& equation) override { + Plog::validate(equation); + } +}; + + +//! Pressure-dependent rate expression where the rate coefficient is expressed +//! as a bivariate Chebyshev polynomial in temperature and pressure. +/*! + * The rate constant can be written as: + * \f[ + * \log k(T,P) = \sum_{t=1}^{N_T} \sum_{p=1}^{N_P} \alpha_{tp} + * \phi_t(\tilde{T}) \phi_p(\tilde{P}) + * \f] + * where \f$\alpha_{tp}\f$ are the constants defining the rate, \f$\phi_n(x)\f$ + * is the Chebyshev polynomial of the first kind of degree *n* evaluated at + * *x*, and + * \f[ + * \tilde{T} \equiv \frac{2T^{-1} - T_\mathrm{min}^{-1} - T_\mathrm{max}^{-1}} + * {T_\mathrm{max}^{-1} - T_\mathrm{min}^{-1}} + * \f] + * \f[ + * \tilde{P} \equiv \frac{2 \log P - \log P_\mathrm{min} - \log P_\mathrm{max}} + * {\log P_\mathrm{max} - \log P_\mathrm{min}} + * \f] + * are reduced temperature and reduced pressures which map the ranges + * \f$ (T_\mathrm{min}, T_\mathrm{max}) \f$ and + * \f$ (P_\mathrm{min}, P_\mathrm{max}) \f$ to (-1, 1). + * + * A Chebyshev rate expression is specified in terms of the coefficient matrix + * \f$ \alpha \f$ and the temperature and pressure ranges. Note that the + * Chebyshev polynomials are not defined outside the interval (-1,1), and + * therefore extrapolation of rates outside the range of temperatures and + * pressures for which they are defined is strongly discouraged. + */ +class ChebyshevRate3 final : public ReactionRate, public Chebyshev +{ +public: + //! Default constructor. + ChebyshevRate3() {} + + //! Constructor directly from coefficient array + /* + * @param Tmin Minimum temperature [K] + * @param Tmax Maximum temperature [K] + * @param Pmin Minimum pressure [Pa] + * @param Pmax Maximum pressure [Pa] + * @param coeffs Coefficient array dimensioned `nT` by `nP` where `nT` and + * `nP` are the number of temperatures and pressures used in the fit, + * respectively. + */ + ChebyshevRate3(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs); + + //! Constructor using AnyMap content + //! @param node AnyMap containing rate information + //! @param rate_units unit definitions used for rate information + ChebyshevRate3(const AnyMap& node, const Units& rate_units); + + //! Constructor using AnyMap content + //! @param node AnyMap containing rate information + ChebyshevRate3(const AnyMap& node); + + virtual std::string type() const override { return "ChebyshevRate"; } + + virtual void setParameters(const AnyMap& node, const Units& rate_units) override; + virtual void getParameters(AnyMap& rateNode, + const Units& rate_units) const override; + + //! Update information specific to reaction + static bool usesUpdate() { return true; } + + virtual void update(const ChebyshevData& shared_data, double concm=0.) override { + update_C(&shared_data.m_log10P); + } + + virtual double eval(const ChebyshevData& shared_data, + double concm=0.) const override { + return updateRC(0., shared_data.m_recipT); + } + + virtual void validate(const std::string& equation) override; }; @@ -182,8 +435,12 @@ class CustomFunc1Rate final : public ReactionRate virtual std::string type() const override { return "custom-function"; } + virtual void setParameters(const AnyMap& node, const Units& rate_units) override { + units = rate_units; + } + //! Update information specific to reaction - static bool uses_update() { return false; } + static bool usesUpdate() { return false; } //! Set custom rate /** @@ -192,7 +449,10 @@ class CustomFunc1Rate final : public ReactionRate */ void setRateFunction(shared_ptr f); - virtual double eval(const CustomFunc1Data& shared_data) const override; + virtual double eval(const CustomFunc1Data& shared_data, + double concm=0.) const override; + + virtual void validate(const std::string& equation) override {} protected: shared_ptr m_ratefunc; diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 5f526ebf220..2edc6dbe7a2 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -1,5 +1,8 @@ /** * @file RxnRates.h + * + * @TODO at least part of the content of this file should be transferred + * to ReactionRate.h once the old XML interface is removed after Cantera 2.6 */ // This file is part of Cantera. See License.txt in the top-level directory or @@ -16,6 +19,7 @@ namespace Cantera class Array2D; class AnyValue; +class AnyMap; class UnitSystem; class Units; class AnyMap; @@ -43,7 +47,14 @@ class Arrhenius /// @param E Activation energy in temperature units. Kelvin. Arrhenius(doublereal A, doublereal b, doublereal E); - //! Run object setup based on AnyMap node information + //! Constructor based on AnyMap content + Arrhenius(const AnyValue& rate, + const UnitSystem& units, const Units& rate_units); + + //! Perform object setup based on AnyMap node information + //! @param node AnyMap containing rate information + //! @param units unit system + //! @param rate_units unit definitions specific to rate information void setParameters(const AnyValue& rate, const UnitSystem& units, const Units& rate_units); @@ -319,11 +330,23 @@ class Plog { public: //! Default constructor. - Plog() {} + Plog(); //! Constructor from Arrhenius rate expressions at a set of pressures explicit Plog(const std::multimap& rates); + //! Perform object setup based on AnyMap node information + //! @param rates vector of AnyMap containing rate information + //! @param units unit system + //! @param rate_units unit definitions specific to rate information + void setParameters(const std::vector& rates, + const UnitSystem& units, const Units& rate_units); + + void getParameters(AnyMap& rateNode, const Units& rate_units) const; + + //! Set up Plog object + void setup(const std::multimap& rates); + //! Update concentration-dependent parts of the rate coefficient. //! @param c natural log of the pressure in Pa void update_C(const doublereal* c) { @@ -439,11 +462,11 @@ class Plog * therefore extrapolation of rates outside the range of temperatures and * pressures for which they are defined is strongly discouraged. */ -class ChebyshevRate +class Chebyshev { public: //! Default constructor. - ChebyshevRate() {} + Chebyshev() : nP_(0), nT_(0) {} //! Constructor directly from coefficient array /* @@ -455,12 +478,24 @@ class ChebyshevRate * `nP` are the number of temperatures and pressures used in the fit, * respectively. */ - ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, + Chebyshev(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs); + + //! Perform object setup based on AnyMap node information + //! @param node AnyMap containing rate information + //! @param units unit system + //! @param rate_units unit definitions specific to rate information + void setParameters(const AnyMap& node, + const UnitSystem& units, const Units& rate_units); + void getParameters(AnyMap& rateNode, const Units& rate_units) const; + + //! Set up Chebyshev object + void setup(double Tmin, double Tmax, double Pmin, double Pmax, const Array2D& coeffs); //! Update concentration-dependent parts of the rate coefficient. //! @param c base-10 logarithm of the pressure in Pa - void update_C(const doublereal* c) { + void update_C(const double* c) { double Pr = (2 * c[0] + PrNum_) * PrDen_; double Cnm1 = Pr; double Cn = 1; @@ -483,7 +518,7 @@ class ChebyshevRate * * This function returns the actual value of the rate constant. */ - doublereal updateRC(doublereal logT, doublereal recipT) const { + double updateRC(double logT, double recipT) const { double Tr = (2 * recipT + TrNum_) * TrDen_; double Cnm1 = Tr; double Cn = 1; @@ -687,6 +722,34 @@ class BMSurfaceArrhenius }; + +//! Pressure-dependent rate expression where the rate coefficient is expressed +//! as a bivariate Chebyshev polynomial in temperature and pressure. +/*! + * @deprecated Renamed to Chebyshev. Behavior will change after Cantera 2.6. + * For future behavior, refer to ChebyshevRate3. + */ +class ChebyshevRate : public Chebyshev +{ +public: + //! Default constructor. + ChebyshevRate(); + + //! Constructor directly from coefficient array + /* + * @param Tmin Minimum temperature [K] + * @param Tmax Maximum temperature [K] + * @param Pmin Minimum pressure [Pa] + * @param Pmax Maximum pressure [Pa] + * @param coeffs Coefficient array dimensioned `nT` by `nP` where `nT` and + * `nP` are the number of temperatures and pressures used in the fit, + * respectively. + */ + ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs); + +}; + } #endif diff --git a/include/cantera/kinetics/reaction_defs.h b/include/cantera/kinetics/reaction_defs.h index 7e3050b6427..ed9697ad297 100644 --- a/include/cantera/kinetics/reaction_defs.h +++ b/include/cantera/kinetics/reaction_defs.h @@ -59,11 +59,6 @@ const int PLOG_RXN = 5; */ const int CHEBYSHEV_RXN = 6; -/** - * A custom Python reaction - */ -const int CUSTOMPY_RXN = 99; - /** * A chemical activation reaction. For these reactions, the rate falls * off as the pressure increases, due to collisional stabilization of diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 028dd7318a0..4e45b527195 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -356,21 +356,44 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxReactionRateBase "Cantera::ReactionRateBase": CxxReactionRateBase() string type() + void update(double) except +translate_exception + void update(double, double) except +translate_exception double eval(double) except +translate_exception double eval(double, double) except +translate_exception double ddT(double) except +translate_exception double ddT(double, double) except +translate_exception - - cdef cppclass CxxCustomFunc1Rate "Cantera::CustomFunc1Rate" (CxxReactionRateBase): - CxxCustomFunc1Rate() - void setRateFunction(shared_ptr[CxxFunc1]) except +translate_exception + CxxAnyMap parameters() except +translate_exception cdef cppclass CxxArrheniusRate "Cantera::ArrheniusRate" (CxxReactionRateBase): CxxArrheniusRate() + CxxArrheniusRate(CxxAnyMap) except +translate_exception CxxArrheniusRate(double, double, double) double preExponentialFactor() double temperatureExponent() - double activationEnergy_R() + double activationEnergy() + cbool allow_negative_pre_exponential_factor + + cdef cppclass CxxPlogRate "Cantera::PlogRate" (CxxReactionRateBase): + CxxPlogRate() + CxxPlogRate(CxxAnyMap) except +translate_exception + CxxPlogRate(multimap[double, CxxArrhenius]) + vector[pair[double, CxxArrhenius]] rates() + + cdef cppclass CxxChebyshevRate3 "Cantera::ChebyshevRate3" (CxxReactionRateBase): + CxxChebyshevRate3() + CxxChebyshevRate3(CxxAnyMap) except +translate_exception + CxxChebyshevRate3(double, double, double, double, CxxArray2D) + double Tmin() + double Tmax() + double Pmin() + double Pmax() + size_t nPressure() + size_t nTemperature() + vector[double]& coeffs() + + cdef cppclass CxxCustomFunc1Rate "Cantera::CustomFunc1Rate" (CxxReactionRateBase): + CxxCustomFunc1Rate() + void setRateFunction(shared_ptr[CxxFunc1]) except +translate_exception cdef cppclass CxxReaction "Cantera::Reaction": CxxReaction() @@ -391,14 +414,10 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cbool duplicate cbool allow_nonreactant_orders cbool allow_negative_orders + cbool usesLegacy() - cdef cppclass CxxReaction2 "Cantera::Reaction2" (CxxReaction): - CxxReaction2() - shared_ptr[CxxReactionRateBase] rate() - void setRate(shared_ptr[CxxReactionRateBase]) - - cdef cppclass CxxElementaryReaction "Cantera::ElementaryReaction" (CxxReaction): - CxxElementaryReaction() + cdef cppclass CxxElementaryReaction2 "Cantera::ElementaryReaction2" (CxxReaction): + CxxElementaryReaction2() CxxArrhenius rate cbool allow_negative_pre_exponential_factor @@ -409,8 +428,8 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": Composition efficiencies double default_efficiency - cdef cppclass CxxThreeBodyReaction "Cantera::ThreeBodyReaction" (CxxElementaryReaction): - CxxThreeBodyReaction() + cdef cppclass CxxThreeBodyReaction2 "Cantera::ThreeBodyReaction2" (CxxElementaryReaction2): + CxxThreeBodyReaction2() CxxThirdBody third_body cdef cppclass CxxFalloff "Cantera::Falloff": @@ -441,11 +460,11 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": void update_C(double*) double updateRC(double, double) - cdef cppclass CxxPlogReaction "Cantera::PlogReaction" (CxxReaction): + cdef cppclass CxxPlogReaction2 "Cantera::PlogReaction2" (CxxReaction): CxxPlog rate - cdef cppclass CxxChebyshevRate "Cantera::ChebyshevRate": - CxxChebyshevRate(double, double, double, double, CxxArray2D) + cdef cppclass CxxChebyshev "Cantera::Chebyshev": + CxxChebyshev(double, double, double, double, CxxArray2D) double Tmin() double Tmax() double Pmin() @@ -456,15 +475,8 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": void update_C(double*) double updateRC(double, double) - cdef cppclass CxxChebyshevReaction "Cantera::ChebyshevReaction" (CxxReaction): - CxxChebyshevRate rate - - cdef cppclass CxxCustomFunc1Reaction "Cantera::CustomFunc1Reaction" (CxxReaction2): - CxxCustomFunc1Reaction() - - cdef cppclass CxxTestReaction "Cantera::TestReaction" (CxxReaction2): - CxxTestReaction() - cbool allow_negative_pre_exponential_factor + cdef cppclass CxxChebyshevReaction2 "Cantera::ChebyshevReaction2" (CxxReaction): + CxxChebyshev rate cdef cppclass CxxBlowersMasel "Cantera::BlowersMasel": CxxBlowersMasel() @@ -487,7 +499,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": double E double m - cdef cppclass CxxInterfaceReaction "Cantera::InterfaceReaction" (CxxElementaryReaction): + cdef cppclass CxxInterfaceReaction "Cantera::InterfaceReaction" (CxxElementaryReaction2): stdmap[string, CxxCoverageDependency] coverage_deps cbool is_sticking_coefficient cbool use_motz_wise_correction @@ -499,6 +511,27 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cbool use_motz_wise_correction string sticking_species + cdef cppclass CxxReaction3 "Cantera::Reaction3" (CxxReaction): + CxxReaction3() + shared_ptr[CxxReactionRateBase] rate() + void setRate(shared_ptr[CxxReactionRateBase]) + + cdef cppclass CxxElementaryReaction3 "Cantera::ElementaryReaction3" (CxxReaction3): + CxxElementaryReaction3() + + cdef cppclass CxxThreeBodyReaction3 "Cantera::ThreeBodyReaction3" (CxxElementaryReaction3): + CxxThreeBodyReaction3() + shared_ptr[CxxThirdBody] thirdBody() + + cdef cppclass CxxPlogReaction3 "Cantera::PlogReaction3" (CxxReaction3): + CxxPlogReaction3() + + cdef cppclass CxxChebyshevReaction3 "Cantera::ChebyshevReaction3" (CxxReaction3): + CxxChebyshevReaction3() + + cdef cppclass CxxCustomFunc1Reaction "Cantera::CustomFunc1Reaction" (CxxReaction3): + CxxCustomFunc1Reaction() + cdef extern from "cantera/kinetics/FalloffFactory.h" namespace "Cantera": cdef shared_ptr[CxxFalloff] CxxNewFalloff "Cantera::newFalloff" (string, vector[double]) except +translate_exception @@ -1117,15 +1150,25 @@ cdef class _ReactionRate: cdef shared_ptr[CxxReactionRateBase] _base cdef CxxReactionRateBase* base -cdef class CustomRate(_ReactionRate): - cdef CxxCustomFunc1Rate* rate - cdef Func1 _rate_func - cdef class ArrheniusRate(_ReactionRate): cdef CxxArrheniusRate* rate @staticmethod cdef wrap(shared_ptr[CxxReactionRateBase]) +cdef class PlogRate(_ReactionRate): + cdef CxxPlogRate* rate + @staticmethod + cdef wrap(shared_ptr[CxxReactionRateBase]) + +cdef class ChebyshevRate(_ReactionRate): + cdef CxxChebyshevRate3* rate + @staticmethod + cdef wrap(shared_ptr[CxxReactionRateBase]) + +cdef class CustomRate(_ReactionRate): + cdef CxxCustomFunc1Rate* rate + cdef Func1 _rate_func + cdef class Reaction: cdef shared_ptr[CxxReaction] _reaction cdef CxxReaction* reaction diff --git a/interfaces/cython/cantera/examples/kinetics/blowers_masel.py b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py index 90db4174c69..de7fff54d94 100644 --- a/interfaces/cython/cantera/examples/kinetics/blowers_masel.py +++ b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py @@ -24,7 +24,7 @@ #Create an elementary reaction O+H2<=>H+OH r1 = ct.ElementaryReaction({'O':1, 'H2':1}, {'H':1, 'OH':1}) -r1.rate = ct.Arrhenius(3.87e1, 2.7, 6260*1000*4.184) +r1.rate = ct.ArrheniusRate(3.87e1, 2.7, 6260*1000*4.184) #Create a Blowers-Masel reaction O+H2<=>H+OH r2 = ct.BlowersMaselReaction({'O':1, 'H2':1}, {'H':1, 'OH':1}) diff --git a/interfaces/cython/cantera/examples/kinetics/custom_reactions.py b/interfaces/cython/cantera/examples/kinetics/custom_reactions.py index 267d746569d..3d39a3dcb9b 100644 --- a/interfaces/cython/cantera/examples/kinetics/custom_reactions.py +++ b/interfaces/cython/cantera/examples/kinetics/custom_reactions.py @@ -29,24 +29,9 @@ gas1 = ct.Solution(thermo='ideal-gas', kinetics='gas', species=species, reactions=custom_reactions) -# construct test reactions: replace all elementary reactions with alternatives -test_reactions = [] -for r in reactions: +# old framework - use xml input - if r.reaction_type == "elementary": - A = r.rate.pre_exponential_factor - b = r.rate.temperature_exponent - Ea = r.rate.activation_energy - r_new = ct.TestReaction(equation=r.equation, - rate={'A': A, 'b': b, 'Ea': Ea}, - kinetics=gas0) - else: - r_new = r - - test_reactions.append(r_new) - -gas2 = ct.Solution(thermo='ideal-gas', kinetics='gas', - species=species, reactions=test_reactions) +gas2 = ct.Solution(mech.replace('.yaml', '.xml')) # construct test case - simulate ignition @@ -75,7 +60,7 @@ def ignition(gas): for i in range(repeat): sim0 += ignition(gas0) sim0 /= repeat -print('- Original mechanism: ' +print('- New framework (YAML): ' '{0:.2f} ms (T_final={1:.2f})'.format(sim0, gas0.T)) sim1 = 0 @@ -90,6 +75,6 @@ def ignition(gas): for i in range(repeat): sim2 += ignition(gas2) sim2 /= repeat -print('- Alternative reactions: ' +print('- Old framework (XML): ' '{0:.2f} ms (T_final={1:.2f}) ... ' '{2:+.2f}%'.format(sim2, gas2.T, 100 * sim2 / sim0 - 100)) diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index 7c31a16319c..cd025c8587e 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -6,6 +6,315 @@ cdef dict _reaction_class_registry = {} +cdef class _ReactionRate: + """ + Base class for ReactionRate objects. + + ReactionRate objects are used to calculate reaction rates and are associated + with a Reaction object. + """ + + def __repr__(self): + return f"<{pystr(self.base.type())} at {id(self):0x}>" + + def __call__(self, double temperature, pressure=None): + """ + Evaluate rate expression based on temperature and pressure. For rate + expressions that are dependent of pressure, an omission of pressure + will raise an exception. + """ + if pressure: + self.base.update(temperature, pressure) + return self.base.eval(temperature, pressure) + else: + self.base.update(temperature) + return self.base.eval(temperature) + + def ddT(self, double temperature, pressure=None): + """ + Evaluate derivative of rate expression with respect to temperature. + """ + if pressure: + self.base.update(temperature, pressure) + return self.base.ddT(temperature, pressure) + else: + self.base.update(temperature) + return self.base.ddT(temperature) + + property input_data: + """ + Get input data for this reaction rate with its current parameter values. + """ + def __get__(self): + return anymap_to_dict(self.base.parameters()) + + +cdef class ArrheniusRate(_ReactionRate): + r""" + A reaction rate coefficient which depends on temperature only and follows + the modified Arrhenius form: + + .. math:: + + k_f = A T^b \exp{-\tfrac{E}{RT}} + + where *A* is the `pre_exponential_factor`, *b* is the `temperature_exponent`, + and *Ea* is the `activation_energy`. + """ + def __cinit__(self, A=None, b=None, Ea=None, input_data=None, init=True): + + if init: + if isinstance(input_data, dict): + self._base.reset(new CxxArrheniusRate(dict_to_anymap(input_data))) + elif all([arg is not None for arg in [A, b, Ea]]): + self._base.reset(new CxxArrheniusRate(A, b, Ea)) + elif all([arg is None for arg in [A, b, Ea, input_data]]): + self._base.reset(new CxxArrheniusRate(dict_to_anymap({}))) + elif input_data: + raise TypeError("Invalid parameter 'input_data'") + else: + raise TypeError("Invalid parameters 'A', 'b' or 'Ea'") + self.base = self._base.get() + self.rate = (self.base) + + @staticmethod + cdef wrap(shared_ptr[CxxReactionRateBase] rate): + """ + Wrap a C++ ReactionRateBase object with a Python object. + """ + # wrap C++ reaction + cdef ArrheniusRate arr + arr = ArrheniusRate(init=False) + arr._base = rate + arr.base = arr._base.get() + arr.rate = (arr.base) + return arr + + property pre_exponential_factor: + """ + The pre-exponential factor *A* in units of m, kmol, and s raised to + powers depending on the reaction order. + """ + def __get__(self): + return self.rate.preExponentialFactor() + + property temperature_exponent: + """ + The temperature exponent *b*. + """ + def __get__(self): + return self.rate.temperatureExponent() + + property activation_energy: + """ + The activation energy *E* [J/kmol]. + """ + def __get__(self): + return self.rate.activationEnergy() + + property allow_negative_pre_exponential_factor: + """ + Get/Set whether the rate coefficient is allowed to have a negative + pre-exponential factor. + """ + def __get__(self): + return self.rate.allow_negative_pre_exponential_factor + def __set__(self, allow): + self.rate.allow_negative_pre_exponential_factor = allow + + +cdef class PlogRate(_ReactionRate): + r""" + A pressure-dependent reaction rate parameterized by logarithmically + interpolating between Arrhenius rate expressions at various pressures. + """ + def __cinit__(self, rates=None, input_data=None, init=True): + + if init and isinstance(rates, list): + self.rates = rates + + elif init: + if isinstance(input_data, dict): + self._base.reset(new CxxPlogRate(dict_to_anymap(input_data))) + elif rates is None: + self._base.reset(new CxxPlogRate(dict_to_anymap({}))) + elif input_data: + raise TypeError("Invalid parameter 'input_data'") + else: + raise TypeError("Invalid parameter 'rates'") + self.base = self._base.get() + self.rate = (self.base) + + @staticmethod + cdef wrap(shared_ptr[CxxReactionRateBase] rate): + """ + Wrap a C++ ReactionRateBase object with a Python object. + """ + # wrap C++ reaction + cdef PlogRate arr + arr = PlogRate(init=False) + arr._base = rate + arr.base = arr._base.get() + arr.rate = (arr.base) + return arr + + property rates: + """ + Get/Set the rate coefficients for this reaction, which are given as a + list of (pressure, `Arrhenius`) tuples. + """ + def __get__(self): + rates = [] + cdef vector[pair[double, CxxArrhenius]] cxxrates = self.rate.rates() + cdef pair[double, CxxArrhenius] p_rate + for p_rate in cxxrates: + rates.append((p_rate.first, copyArrhenius(&p_rate.second))) + return rates + + def __set__(self, rates): + cdef multimap[double, CxxArrhenius] ratemap + cdef Arrhenius rate + cdef pair[double, CxxArrhenius] item + for p, rate in rates: + item.first = p + item.second = deref(rate.rate) + ratemap.insert(item) + + self._base.reset(new CxxPlogRate(ratemap)) + self.base = self._base.get() + self.rate = (self.base) + + +cdef class ChebyshevRate(_ReactionRate): + r""" + A pressure-dependent reaction rate parameterized by a bivariate Chebyshev + polynomial in temperature and pressure. + """ + def __cinit__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, data=None, + input_data=None, init=True): + + if init: + if isinstance(input_data, dict): + self._base.reset(new CxxChebyshevRate3(dict_to_anymap(input_data))) + elif all([arg is not None for arg in [Tmin, Tmax, Pmin, Pmax, data]]): + self._setup(Tmin, Tmax, Pmin, Pmax, data) + return + elif all([arg is None + for arg in [Tmin, Tmax, Pmin, Pmax, data, input_data]]): + self._base.reset(new CxxChebyshevRate3(dict_to_anymap({}))) + elif input_data: + raise TypeError("Invalid parameter 'input_data'") + else: + raise TypeError("Invalid parameters") + self.base = self._base.get() + self.rate = (self.base) + + def _setup(self, Tmin, Tmax, Pmin, Pmax, coeffs): + """ + Simultaneously set values for `Tmin`, `Tmax`, `Pmin`, `Pmax`, and + `coeffs`. + """ + cdef CxxArray2D data + data.resize(len(coeffs), len(coeffs[0])) + cdef double value + cdef int i + cdef int j + for i,row in enumerate(coeffs): + for j,value in enumerate(row): + CxxArray2D_set(data, i, j, value) + + self._base.reset(new CxxChebyshevRate3(Tmin, Tmax, Pmin, Pmax, data)) + self.base = self._base.get() + self.rate = (self.base) + + @staticmethod + cdef wrap(shared_ptr[CxxReactionRateBase] rate): + """ + Wrap a C++ ReactionRateBase object with a Python object. + """ + # wrap C++ reaction + cdef ChebyshevRate arr + arr = ChebyshevRate(init=False) + arr._base = rate + arr.base = arr._base.get() + arr.rate = (arr.base) + return arr + + property Tmin: + """ Minimum temperature [K] for the Chebyshev fit """ + def __get__(self): + return self.rate.Tmin() + + property Tmax: + """ Maximum temperature [K] for the Chebyshev fit """ + def __get__(self): + return self.rate.Tmax() + + property Pmin: + """ Minimum pressure [Pa] for the Chebyshev fit """ + def __get__(self): + return self.rate.Pmin() + + property Pmax: + """ Maximum pressure [Pa] for the Chebyshev fit """ + def __get__(self): + return self.rate.Pmax() + + property n_pressure: + """ Number of pressures over which the Chebyshev fit is computed """ + def __get__(self): + return self.rate.nPressure() + + property n_temperature: + """ Number of temperatures over which the Chebyshev fit is computed """ + def __get__(self): + return self.rate.nTemperature() + + property coeffs: + """ + 2D array of Chebyshev coefficients of size `(n_temperature, n_pressure)`. + """ + def __get__(self): + c = np.fromiter(self.rate.coeffs(), np.double) + return c.reshape((self.rate.nTemperature(), self.rate.nPressure())) + + +cdef class CustomRate(_ReactionRate): + r""" + A custom rate coefficient which depends on temperature only. + + The simplest way to create a `CustomRate` object is to use a lambda function, + for example:: + + rr = CustomRate(lambda T: 38.7 * T**2.7 * exp(-3150.15/T)) + """ + def __cinit__(self, k=None, init=True): + + if init: + self._base.reset(new CxxCustomFunc1Rate()) + self.base = self._base.get() + self.rate = (self.base) + self.set_rate_function(k) + + def set_rate_function(self, k): + r""" + Set the function describing a custom reaction rate:: + + rr = CustomRate() + rr.set_rate_function(lambda T: 38.7 * T**2.7 * exp(-3150.15/T)) + """ + if k is None: + self._rate_func = None + return + + if isinstance(k, Func1): + self._rate_func = k + else: + self._rate_func = Func1(k) + + self.rate.setRateFunction(self._rate_func._func) + + cdef class Reaction: """ A class which stores data about a reaction and its rate parameterization so @@ -22,9 +331,9 @@ cdef class Reaction: will produce a list of the 325 reactions which make up the GRI 3.0 mechanism:: - R = ct.Reaction.listFromFile('gri30.yaml', gas) - R = ct.Reaction.listFromCti(open('path/to/gri30.cti').read()) - R = ct.Reaction.listFromXml(open('path/to/gri30.xml').read()) + R = ct.Reaction.listFromFile("gri30.yaml", gas) + R = ct.Reaction.listFromCti(open("path/to/gri30.cti").read()) + R = ct.Reaction.listFromXml(open("path/to/gri30.xml").read()) where `gas` is a `Solution` object with the appropriate thermodynamic model, which is the `ideal-gas` model in this case. @@ -56,12 +365,17 @@ cdef class Reaction: R = ct.Reaction.fromCti('''reaction('O + H2 <=> H + OH', [(3.87e4, 'cm3/mol/s'), 2.7, (6260, 'cal/mol')])''') """ - reaction_type = "" + _reaction_type = "" + _has_legacy = False + _hybrid = False # indicate whether legacy implementations are separate or merged - def __cinit__(self, reactants='', products='', init=True, **kwargs): + def __cinit__(self, reactants="", products="", init=True, legacy=False, **kwargs): if init: - self._reaction = CxxNewReaction(stringify((self.reaction_type))) + rxn_type = self._reaction_type + if (not self._hybrid and self._has_legacy) or (self._hybrid and legacy): + rxn_type += "-legacy" + self._reaction = CxxNewReaction(stringify((rxn_type))) self.reaction = self._reaction.get() if reactants: self.reactants = reactants @@ -77,7 +391,13 @@ cdef class Reaction: if not(_reaction_class_registry): def register_subclasses(cls): for c in cls.__subclasses__(): - _reaction_class_registry[getattr(c, 'reaction_type')] = c + rxn_type = getattr(c, "_reaction_type") + if getattr(c, "_hybrid"): + # registry needs to contain both updated and "-legacy" variants + _reaction_class_registry[rxn_type] = c + if getattr(c, "_has_legacy", False): + rxn_type += "-legacy" + _reaction_class_registry[rxn_type] = c register_subclasses(c) # update global reaction class registry @@ -119,19 +439,68 @@ cdef class Reaction: cxx_reaction = CxxNewReaction(deref(CxxGetXmlFromString(stringify(text)))) return Reaction.wrap(cxx_reaction) - @staticmethod - def fromYaml(text, Kinetics kinetics): + @classmethod + def from_dict(cls, data, Kinetics kinetics=None): + """ + Create a `Reaction` object from a dictionary corresponding to its YAML + representation. + + An example for the creation of a Reaction from a dictionary is:: + + rxn = Reaction.from_dict( + {"equation": "O + H2 <=> H + OH", + "rate-constant": {"A": 38.7, "b": 2.7, "Ea": 26191840.0}}, + kinetics=gas) + + In the example, *gas* is a Kinetics (or Solution) object. + + :param data: + A dictionary corresponding to the YAML representation. + :param kinetics: + A `Kinetics` object whose associated phase(s) contain the species + involved in the reaction. + """ + if cls._reaction_type != "": + raise TypeError( + "Class method 'from_dict' was invoked from '{}' but should " + "be called from base class 'Reaction'".format(cls.__name__)) + if kinetics is None: + raise ValueError("A Kinetics object is required.") + + cdef CxxAnyMap any_map = dict_to_anymap(data) + cxx_reaction = CxxNewReaction(any_map, deref(kinetics.kinetics)) + return Reaction.wrap(cxx_reaction) + + @classmethod + def fromYaml(cls, text, Kinetics kinetics=None): """ Create a `Reaction` object from its YAML string representation. + An example for the creation of a Reaction from a YAML string is:: + + rxn = Reaction.fromYaml(''' + equation: O + H2 <=> H + OH + rate-constant: {A: 38.7, b: 2.7, Ea: 6260.0 cal/mol} + ''', kinetics=gas) + + In the example, *gas* is a Kinetics (or Solution) object. + :param text: - The YAML reaction string + The YAML reaction string. :param kinetics: A `Kinetics` object whose associated phase(s) contain the species involved in the reaction. """ - cxx_reaction = CxxNewReaction(AnyMapFromYamlString(stringify(text)), - deref(kinetics.kinetics)) + if cls._reaction_type != "": + raise TypeError( + "Class method 'fromYaml' was invoked from '{}' but should " + "be called from base class 'Reaction'".format(cls.__name__)) + if kinetics is None: + raise ValueError("A Kinetics object is required.") + + cdef CxxAnyMap any_map + any_map = AnyMapFromYamlString(stringify(text)) + cxx_reaction = CxxNewReaction(any_map, deref(kinetics.kinetics)) return Reaction.wrap(cxx_reaction) @staticmethod @@ -280,6 +649,13 @@ cdef class Reaction: def __set__(self, ID): self.reaction.id = stringify(ID) + property reaction_type: + """ + Retrieve the native type name of the reaction. + """ + def __get__(self): + return pystr(self.reaction.type()) + property reversible: """ Get/Set a flag which is `True` if this reaction is reversible or `False` @@ -346,11 +722,21 @@ cdef class Reaction: self.reaction.input.clear() def __repr__(self): - return '<{}: {}>'.format(self.__class__.__name__, self.equation) + return f"<{self.__class__.__name__}: {self.equation}>" def __str__(self): return self.equation + def _deprecation_warning(self, attr, what="property"): + return ("\n{} '{}' to be removed after Cantera 2.6.\nThis {} is moved to " + "the {} object accessed via the 'rate' property." + ).format(what.capitalize(), attr, what, type(self.rate).__name__) + + property uses_legacy: + """Indicate whether reaction uses a legacy implementation""" + def __get__(self): + return self.reaction.usesLegacy() + cdef class Arrhenius: r""" @@ -415,193 +801,199 @@ cdef wrapArrhenius(CxxArrhenius* rate, Reaction reaction): r.reaction = reaction return r -cdef copyArrhenius(CxxArrhenius* rate): - r = Arrhenius(rate.preExponentialFactor(), rate.temperatureExponent(), - rate.activationEnergy_R() * gas_constant) - return r - - -cdef class _ReactionRate: - - def __repr__(self): - return "<{} at {:0x}>".format(pystr(self.base.type()), id(self)) - - def __call__(self, double temperature, pressure=None): - if pressure: - return self.base.eval(temperature, pressure) - else: - return self.base.eval(temperature) - - def ddT(self, double temperature, pressure=None): - if pressure: - return self.base.ddT(temperature, pressure) - else: - return self.base.ddT(temperature) - - -cdef class CustomRate(_ReactionRate): - r""" - A custom rate coefficient which depends on temperature only. - - The simplest way to create a `CustomRate` object is to use a lambda function, - for example:: - - rr = CustomRate(lambda T: 38.7 * T**2.7 * exp(-3150.15/T)) - - Warning: this class is an experimental part of the Cantera API and - may be changed or removed without notice. - """ - def __cinit__(self, k=None, init=True): - - if init: - self._base.reset(new CxxCustomFunc1Rate()) - self.base = self._base.get() - self.rate = (self.base) - self.set_rate_function(k) - - def set_rate_function(self, k): - r""" - Set the function describing a custom reaction rate:: - - rr = CustomRate() - rr.set_rate_function(lambda T: 38.7 * T**2.7 * exp(-3150.15/T)) - """ - if k is None: - self._rate_func = None - return - - if isinstance(k, Func1): - self._rate_func = k - else: - self._rate_func = Func1(k) - - self.rate.setRateFunction(self._rate_func._func) - - -cdef class ArrheniusRate(_ReactionRate): - r""" - A reaction rate coefficient which depends on temperature only and follows - the modified Arrhenius form: - - .. math:: - - k_f = A T^b \exp{-\tfrac{E}{RT}} - - where *A* is the `pre_exponential_factor`, *b* is the `temperature_exponent`, - and *E* is the `activation_energy`. - - Warning: this class is an experimental part of the Cantera API and - may be changed or removed without notice. - """ - def __cinit__(self, A=0, b=0, E=0, init=True): - - if init: - self._base.reset(new CxxArrheniusRate(A, b, E / gas_constant)) - self.base = self._base.get() - self.rate = (self.base) - - @staticmethod - cdef wrap(shared_ptr[CxxReactionRateBase] rate): - """ - Wrap a C++ ReactionRateBase object with a Python object. - """ - # wrap C++ reaction - cdef ArrheniusRate arr - arr = ArrheniusRate(init=False) - arr._base = rate - arr.base = arr._base.get() - arr.rate = (arr.base) - return arr - - property pre_exponential_factor: - """ - The pre-exponential factor *A* in units of m, kmol, and s raised to - powers depending on the reaction order. - """ - def __get__(self): - return self.rate.preExponentialFactor() - - property temperature_exponent: - """ - The temperature exponent *b*. - """ - def __get__(self): - return self.rate.temperatureExponent() - - property activation_energy: - """ - The activation energy *E* [J/kmol]. - """ - def __get__(self): - return self.rate.activationEnergy_R() * gas_constant - +cdef copyArrhenius(CxxArrhenius* rate): + r = Arrhenius(rate.preExponentialFactor(), rate.temperatureExponent(), + rate.activationEnergy_R() * gas_constant) + return r + cdef class ElementaryReaction(Reaction): """ A reaction which follows mass-action kinetics with a modified Arrhenius reaction rate. - An example for the definition of a `ElementaryReaction` object is given as:: + An example for the definition of an `ElementaryReaction` object is given as:: + + rxn = ElementaryReaction( + equation="O + H2 <=> H + OH", + rate={"A": 38.7, "b": 2.7, "Ea": 2.619184e+07}, + kinetics=gas) - rxn = ElementaryReaction(equation='H2 + O <=> H + OH', - rate={'A': 38.7, 'b': 2.7, 'Ea': 2.619184e+07}, - kinetics=gas) + The YAML description corresponding to this reaction is:: + + equation: O + H2 <=> H + OH + rate-constant: {A: 3.87e+04 cm^3/mol/s, b: 2.7, Ea: 6260.0 cal/mol} """ - reaction_type = "elementary" + _reaction_type = "elementary" + _has_legacy = True + _hybrid = True + + cdef CxxElementaryReaction3* er(self): + if self.uses_legacy: + raise AttributeError("Incorrect accessor for updated implementation") + return self.reaction + + cdef CxxElementaryReaction2* er2(self): + if not self.uses_legacy: + raise AttributeError("Incorrect accessor for legacy implementation") + return self.reaction def __init__(self, equation=None, rate=None, Kinetics kinetics=None, - init=True, **kwargs): + init=True, legacy=False, **kwargs): if init and equation and kinetics: + rxn_type = self._reaction_type + if legacy: + rxn_type += "-legacy" + spec = {"equation": equation, "type": rxn_type} if isinstance(rate, dict): - coeffs = ['{}: {}'.format(k, v) for k, v in rate.items()] - elif isinstance(rate, Arrhenius) or rate is None: - coeffs = ['{}: 0.'.format(k) for k in ['A', 'b', 'Ea']] + spec["rate-constant"] = rate + elif legacy and (isinstance(rate, Arrhenius) or rate is None): + spec["rate-constant"] = dict.fromkeys(["A", "b", "Ea"], 0.) + elif rate is None: + pass + elif not legacy and isinstance(rate, (Arrhenius, ArrheniusRate)): + pass else: raise TypeError("Invalid rate definition") - rate_def = '{{{}}}'.format(', '.join(coeffs)) - yaml = '{{equation: {}, rate-constant: {}, type: {}}}'.format( - equation, rate_def, self.reaction_type) - self._reaction = CxxNewReaction(AnyMapFromYamlString(stringify(yaml)), + self._reaction = CxxNewReaction(dict_to_anymap(spec), deref(kinetics.kinetics)) self.reaction = self._reaction.get() - if isinstance(rate, Arrhenius): + if legacy and isinstance(rate, Arrhenius): self.rate = rate + elif not legacy and isinstance(rate, (ArrheniusRate, Arrhenius)): + self.rate = rate + + cdef _legacy_set_rate(self, Arrhenius rate): + cdef CxxElementaryReaction2* r = self.er2() + r.rate = deref(rate.rate) property rate: - """ Get/Set the `Arrhenius` rate coefficient for this reaction. """ + """ Get/Set the `ArrheniusRate` rate coefficient for this reaction. """ def __get__(self): - cdef CxxElementaryReaction* r = self.reaction - return wrapArrhenius(&(r.rate), self) - def __set__(self, Arrhenius rate): - cdef CxxElementaryReaction* r = self.reaction - r.rate = deref(rate.rate) + if self.uses_legacy: + return wrapArrhenius(&(self.er2().rate), self) + + return ArrheniusRate.wrap(self.er().rate()) + def __set__(self, rate): + if self.uses_legacy: + self._legacy_set_rate(rate) + return + + cdef ArrheniusRate rate3 + if isinstance(rate, ArrheniusRate): + rate3 = rate + elif isinstance(rate, Arrhenius): + warnings.warn("Setting the rate using an 'Arrhenius' object is " + "deprecated and will be removed after Cantera 2.6. The argument " + "type is replaceable by 'ArrheniusRate'.", DeprecationWarning) + rate3 = ArrheniusRate(rate.pre_exponential_factor, + rate.temperature_exponent, + rate.activation_energy) + else: + raise TypeError("Invalid rate definition") + self.er().setRate(rate3._base) property allow_negative_pre_exponential_factor: """ Get/Set whether the rate coefficient is allowed to have a negative pre-exponential factor. + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `ArrheniusRate.allow_negative_pre_exponential_factor`. """ def __get__(self): - cdef CxxElementaryReaction* r = self.reaction - return r.allow_negative_pre_exponential_factor + if self.uses_legacy: + return self.er2().allow_negative_pre_exponential_factor + + attr = "allow_negative_pre_exponential_factor" + warnings.warn(self._deprecation_warning(attr), DeprecationWarning) + return self.rate.allow_negative_pre_exponential_factor def __set__(self, allow): - cdef CxxElementaryReaction* r = self.reaction - r.allow_negative_pre_exponential_factor = allow + if self.uses_legacy: + self.er2().allow_negative_pre_exponential_factor = allow + return + + attr = "allow_negative_pre_exponential_factor" + warnings.warn(self._deprecation_warning(attr), DeprecationWarning) + self.rate.allow_negative_pre_exponential_factor = allow cdef class ThreeBodyReaction(ElementaryReaction): """ A reaction with a non-reacting third body "M" that acts to add or remove energy from the reacting species. + + An example for the definition of an `ThreeBodyReaction` object is given as:: + + rxn = ThreeBodyReaction( + equation="2 O + M <=> O2 + M", + type="three-body", + rate={"A": 1.2e+17, "b": -1.0, "Ea": 0.0}, + efficiencies={"H2": 2.4, "H2O": 15.4, "AR": 0.83}, + kinetics=gas) + + The YAML description corresponding to this reaction is:: + + equation: 2 O + M <=> O2 + M + 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" + _reaction_type = "three-body" + _has_legacy = True + _hybrid = True + + cdef CxxThreeBodyReaction3* tbr(self): + if self.uses_legacy: + raise AttributeError("Incorrect accessor for updated implementation") + return self.reaction + + cdef CxxThreeBodyReaction2* tbr2(self): + if not self.uses_legacy: + raise AttributeError("Incorrect accessor for legacy implementation") + return self.reaction + + cdef CxxThirdBody* thirdbody(self): + if self.uses_legacy: + return &(self.tbr2().third_body) + return (self.tbr().thirdBody().get()) - cdef CxxThreeBodyReaction* tbr(self): - return self.reaction + def __init__(self, equation=None, rate=None, efficiencies=None, + Kinetics kinetics=None, legacy=False, init=True, **kwargs): + + if init and equation and kinetics: + + rxn_type = self._reaction_type + if legacy: + rxn_type += "-legacy" + spec = {"equation": equation, "type": rxn_type} + if isinstance(rate, dict): + spec["rate-constant"] = rate + elif legacy and (isinstance(rate, Arrhenius) or rate is None): + spec["rate-constant"] = dict.fromkeys(["A", "b", "Ea"], 0.) + elif rate is None: + pass + elif not legacy and isinstance(rate, (Arrhenius, 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() + + if legacy and isinstance(rate, Arrhenius): + self.rate = rate + elif not legacy and isinstance(rate, (Arrhenius, ArrheniusRate)): + self.rate = rate property efficiencies: """ @@ -610,9 +1002,9 @@ cdef class ThreeBodyReaction(ElementaryReaction): efficiencies. """ def __get__(self): - return comp_map_to_dict(self.tbr().third_body.efficiencies) + return comp_map_to_dict(self.thirdbody().efficiencies) def __set__(self, eff): - self.tbr().third_body.efficiencies = comp_map(eff) + self.thirdbody().efficiencies = comp_map(eff) property default_efficiency: """ @@ -620,16 +1012,16 @@ cdef class ThreeBodyReaction(ElementaryReaction): species used for species not in `efficiencies`. """ def __get__(self): - return self.tbr().third_body.default_efficiency + return self.thirdbody().default_efficiency def __set__(self, default_eff): - self.tbr().third_body.default_efficiency = default_eff + self.thirdbody().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.tbr().third_body.efficiency(stringify(species)) + return self.thirdbody().efficiency(stringify(species)) cdef class Falloff: @@ -727,7 +1119,7 @@ cdef class FalloffReaction(Reaction): A reaction that is first-order in [M] at low pressure, like a third-body reaction, but zeroth-order in [M] as pressure increases. """ - reaction_type = "falloff" + _reaction_type = "falloff" cdef CxxFalloffReaction* frxn(self): return self.reaction @@ -804,44 +1196,143 @@ cdef class ChemicallyActivatedReaction(FalloffReaction): that the forward rate constant is written as being proportional to the low- pressure rate constant. """ - reaction_type = "chemically-activated" + _reaction_type = "chemically-activated" cdef class PlogReaction(Reaction): """ A pressure-dependent reaction parameterized by logarithmically interpolating between Arrhenius rate expressions at various pressures. + + An example for the definition of a `PlogReaction` object is given as:: + + rxn = PlogReaction( + equation="H2 + O2 <=> 2 OH", + rate=[(1013.25, Arrhenius(1.2124e+16, -0.5779, 45491376.8)), + (101325., Arrhenius(4.9108e+31, -4.8507, 103649395.2)), + (1013250., Arrhenius(1.2866e+47, -9.0246, 166508556.0)), + (10132500., Arrhenius(5.9632e+56, -11.529, 220076726.4))], + kinetics=gas) + + The YAML description corresponding to this reaction is:: + + equation: H2 + O2 <=> 2 OH + type: pressure-dependent-Arrhenius + rate-constants: + - {P: 0.01 atm, A: 1.2124e+16, b: -0.5779, Ea: 1.08727e+04 cal/mol} + - {P: 1.0 atm, A: 4.9108e+31, b: -4.8507, Ea: 2.47728e+04 cal/mol} + - {P: 10.0 atm, A: 1.2866e+47, b: -9.0246, Ea: 3.97965e+04 cal/mol} + - {P: 100.0 atm, A: 5.9632e+56, b: -11.529, Ea: 5.25996e+04 cal/mol}or. """ - reaction_type = "pressure-dependent-Arrhenius" + _reaction_type = "pressure-dependent-Arrhenius" + _has_legacy = True + _hybrid = True + + cdef CxxPlogReaction3* pr(self): + if self.uses_legacy: + raise AttributeError("Incorrect accessor for updated implementation") + return self.reaction + + cdef CxxPlogReaction2* cp2(self): + if not self.uses_legacy: + raise AttributeError("Incorrect accessor for legacy implementation") + return self.reaction + + def __init__(self, equation=None, rate=None, Kinetics kinetics=None, + init=True, legacy=False, **kwargs): + + if init and equation and kinetics: + + rxn_type = self._reaction_type + if legacy: + rxn_type += "-legacy" + spec = {"equation": equation, "type": rxn_type} + if legacy and isinstance(rate, list): + rates = [] + for r in rate: + rates.append({ + "P": r[0], + "A": r[1].pre_exponential_factor, + "b": r[1].temperature_exponent, + "Ea": r[1].activation_energy, + }) + spec.update({'rate-constants': rates}) + elif not legacy and (isinstance(rate, (PlogRate, list)) or rate is None): + pass + else: + raise TypeError("Invalid rate definition") + + self._reaction = CxxNewReaction(dict_to_anymap(spec), + deref(kinetics.kinetics)) + self.reaction = self._reaction.get() + + if not legacy and isinstance(rate, PlogRate): + self.rate = rate + elif not legacy and isinstance(rate, list): + self.rate = PlogRate(rate) + + property rate: + """ Get/Set the `PlogRate` rate coefficients for this reaction. """ + def __get__(self): + if self.uses_legacy: + raise AttributeError("Legacy implementation does not use rate property.") + return PlogRate.wrap(self.pr().rate()) + def __set__(self, PlogRate rate): + if self.uses_legacy: + raise AttributeError("Legacy implementation does not use rate property.") + self.pr().setRate(rate._base) + + cdef list _legacy_get_rates(self): + cdef CxxPlogReaction2* r = self.cp2() + cdef vector[pair[double,CxxArrhenius]] cxxrates = r.rate.rates() + cdef pair[double,CxxArrhenius] p_rate + rates = [] + for p_rate in cxxrates: + rates.append((p_rate.first,copyArrhenius(&p_rate.second))) + return rates + + cdef _legacy_set_rates(self, list rates): + cdef multimap[double,CxxArrhenius] ratemap + cdef Arrhenius rate + cdef pair[double,CxxArrhenius] item + for p, rate in rates: + item.first = p + item.second = deref(rate.rate) + ratemap.insert(item) + + cdef CxxPlogReaction2* r = self.cp2() + r.rate = CxxPlog(ratemap) property rates: """ Get/Set the rate coefficients for this reaction, which are given as a list of (pressure, `Arrhenius`) tuples. + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `PlogRate.rates`. """ def __get__(self): - cdef CxxPlogReaction* r = self.reaction - rates = [] - cdef vector[pair[double,CxxArrhenius]] cxxrates = r.rate.rates() - cdef pair[double,CxxArrhenius] p_rate - for p_rate in cxxrates: - rates.append((p_rate.first,copyArrhenius(&p_rate.second))) - return rates - - def __set__(self, rates): - cdef multimap[double,CxxArrhenius] ratemap - cdef Arrhenius rate - cdef pair[double,CxxArrhenius] item - for p,rate in rates: - item.first = p - item.second = deref(rate.rate) - ratemap.insert(item) + if self.uses_legacy: + return self._legacy_get_rates() - cdef CxxPlogReaction* r = self.reaction - r.rate = CxxPlog(ratemap) + warnings.warn(self._deprecation_warning("rates"), DeprecationWarning) + return self.rate.rates - def __call__(self, float T, float P): - cdef CxxPlogReaction* r = self.reaction + def __set__(self, rates): + if self.uses_legacy: + self._legacy_set_rates(rates) + return + + warnings.warn("Property 'rates' to be removed after Cantera 2.6. " + "Setter is replaceable by assigning a new 'PlogRate' object created " + "from rates to the rate property.", DeprecationWarning) + rate_ = self.rate + rate_.rates = rates + self.rate = rate_ + + cdef _legacy_call(self, float T, float P): + cdef CxxPlogReaction2* r = self.cp2() cdef double logT = np.log(T) cdef double recipT = 1/T cdef double logP = np.log(P) @@ -849,65 +1340,209 @@ cdef class PlogReaction(Reaction): r.rate.update_C(&logP) return r.rate.updateRC(logT, recipT) + def __call__(self, float T, float P): + """ + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaceable by call to `rate` property. + """ + if self.uses_legacy: + return self._legacy_call(T, P) + + warnings.warn( + self._deprecation_warning("__call__", "method"), DeprecationWarning) + return self.rate(T, P) + cdef class ChebyshevReaction(Reaction): """ A pressure-dependent reaction parameterized by a bivariate Chebyshev polynomial in temperature and pressure. + + An example for the definition of a `ChebyshevReaction` object is given as:: + + rxn = ChebyshevReaction( + equation="HO2 <=> OH + O", + rate={"Tmin": 290.0, "Tmax": 3000.0, + "Pmin": 1e3, "Pmax": 1e8, + "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]]}, + kinetics=gas) + + The YAML description corresponding to this reaction is:: + + equation: HO2 <=> OH + O + type: Chebyshev + temperature-range: [290.0, 3000.0] + pressure-range: [1.e-03 bar, 10. bar] + 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] """ - reaction_type = "Chebyshev" + _reaction_type = "Chebyshev" + _has_legacy = True + _hybrid = True + + cdef CxxChebyshevReaction3* cr(self): + if self.uses_legacy: + raise AttributeError("Incorrect accessor for updated implementation") + return self.reaction + + cdef CxxChebyshevReaction2* cr2(self): + if not self.uses_legacy: + raise AttributeError("Incorrect accessor for legacy implementation") + return self.reaction + + def __init__(self, equation=None, rate=None, Kinetics kinetics=None, + init=True, legacy=False, **kwargs): + + if init and equation and kinetics: + + rxn_type = self._reaction_type + if legacy: + rxn_type += "-legacy" + spec = {"equation": equation, "type": rxn_type} + if isinstance(rate, dict): + spec["temperature-range"] = [rate["Tmin"], rate["Tmax"]] + spec["pressure-range"] = [rate["Pmin"], rate["Pmax"]] + spec["data"] = rate["data"] + elif not legacy and (isinstance(rate, ChebyshevRate) or rate is None): + pass + else: + raise TypeError("Invalid rate definition") + + self._reaction = CxxNewReaction(dict_to_anymap(spec), + deref(kinetics.kinetics)) + self.reaction = self._reaction.get() + + if not legacy and isinstance(rate, ChebyshevRate): + self.rate = rate + + property rate: + """ Get/Set the `ChebyshevRate` rate coefficients for this reaction. """ + def __get__(self): + if self.uses_legacy: + raise AttributeError("Legacy implementation does not use rate property.") + return ChebyshevRate.wrap(self.cr().rate()) + def __set__(self, ChebyshevRate rate): + if self.uses_legacy: + raise AttributeError("Legacy implementation does not use rate property.") + self.cr().setRate(rate._base) property Tmin: - """ Minimum temperature [K] for the Chebyshev fit """ + """ + Minimum temperature [K] for the Chebyshev fit + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `ChebyshevRate.Tmin`. + """ def __get__(self): - cdef CxxChebyshevReaction* r = self.reaction - return r.rate.Tmin() + if self.uses_legacy: + return self.cr2().rate.Tmin() + + warnings.warn(self._deprecation_warning("Tmin"), DeprecationWarning) + return self.rate.Tmin property Tmax: - """ Maximum temperature [K] for the Chebyshev fit """ + """ + Maximum temperature [K] for the Chebyshev fit + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `ChebyshevRate.Tmax`. + """ def __get__(self): - cdef CxxChebyshevReaction* r = self.reaction - return r.rate.Tmax() + if self.uses_legacy: + return self.cr2().rate.Tmax() + + warnings.warn(self._deprecation_warning("Tmax"), DeprecationWarning) + return self.rate.Tmax property Pmin: - """ Minimum pressure [Pa] for the Chebyshev fit """ + """ + Minimum pressure [Pa] for the Chebyshev fit + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `ChebyshevRate.Pmin`. + """ def __get__(self): - cdef CxxChebyshevReaction* r = self.reaction - return r.rate.Pmin() + if self.uses_legacy: + return self.cr2().rate.Pmin() + + warnings.warn(self._deprecation_warning("Pmin"), DeprecationWarning) + return self.rate.Pmin property Pmax: - """ Maximum pressure [Pa] for the Chebyshev fit """ + """ Maximum pressure [Pa] for the Chebyshev fit + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `ChebyshevRate.Pmax`. + """ def __get__(self): - cdef CxxChebyshevReaction* r = self.reaction - return r.rate.Pmax() + if self.uses_legacy: + return self.cr2().rate.Pmax() + + warnings.warn(self._deprecation_warning("Pmax"), DeprecationWarning) + return self.rate.Pmax property nPressure: - """ Number of pressures over which the Chebyshev fit is computed """ + """ + Number of pressures over which the Chebyshev fit is computed + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `ChebyshevRate.nPressure`. + """ def __get__(self): - cdef CxxChebyshevReaction* r = self.reaction - return r.rate.nPressure() + if self.uses_legacy: + return self.cr2().rate.nPressure() + + warnings.warn(self._deprecation_warning("nPressure"), DeprecationWarning) + return self.rate.n_pressure property nTemperature: - """ Number of temperatures over which the Chebyshev fit is computed """ + """ + Number of temperatures over which the Chebyshev fit is computed + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `ChebyshevRate.nTemperature`. + """ def __get__(self): - cdef CxxChebyshevReaction* r = self.reaction - return r.rate.nTemperature() + if self.uses_legacy: + return self.cr2().rate.nTemperature() + + warnings.warn( + self._deprecation_warning("nTemperature"), DeprecationWarning) + return self.rate.n_temperature + + cdef _legacy_get_coeffs(self): + cdef CxxChebyshevReaction2* r = self.cr2() + c = np.fromiter(r.rate.coeffs(), np.double) + return c.reshape((r.rate.nTemperature(), r.rate.nPressure())) property coeffs: """ - 2D array of Chebyshev coefficients of size `(nTemperature, nPressure)`. + 2D array of Chebyshev coefficients of size `(n_temperature, n_pressure)`. + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by property `ChebyshevRate.coeffs`. """ def __get__(self): - cdef CxxChebyshevReaction* r = self.reaction - c = np.fromiter(r.rate.coeffs(), np.double) - return c.reshape((r.rate.nTemperature(), r.rate.nPressure())) + if self.uses_legacy: + return self._legacy_get_coeffs() - def set_parameters(self, Tmin, Tmax, Pmin, Pmax, coeffs): - """ - Simultaneously set values for `Tmin`, `Tmax`, `Pmin`, `Pmax`, and - `coeffs`. - """ - cdef CxxChebyshevReaction* r = self.reaction + warnings.warn(self._deprecation_warning("coeffs"), DeprecationWarning) + return self.rate.coeffs + + cdef _legacy_set_parameters(self, Tmin, Tmax, Pmin, Pmax, coeffs): + cdef CxxChebyshevReaction2* r = self.cr2() cdef CxxArray2D data data.resize(len(coeffs), len(coeffs[0])) @@ -915,13 +1550,30 @@ cdef class ChebyshevReaction(Reaction): cdef int i cdef int j for i,row in enumerate(coeffs): - for j,value in enumerate(row): + for j, value in enumerate(row): CxxArray2D_set(data, i, j, value) - r.rate = CxxChebyshevRate(Tmin, Tmax, Pmin, Pmax, data) + r.rate = CxxChebyshev(Tmin, Tmax, Pmin, Pmax, data) - def __call__(self, float T, float P): - cdef CxxChebyshevReaction* r = self.reaction + def set_parameters(self, Tmin, Tmax, Pmin, Pmax, coeffs): + """ + Simultaneously set values for `Tmin`, `Tmax`, `Pmin`, `Pmax`, and + `coeffs`. + + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaced by `ChebyshevRate` constructor. + """ + if self.uses_legacy: + return self._legacy_set_parameters(Tmin, Tmax, Pmin, Pmax, coeffs) + + warnings.warn("Method 'set_parameters' to be removed after Cantera 2.6. " + "Method is replaceable by assigning a new 'ChebyshevRate' object to the " + "rate property.", DeprecationWarning) + self.rate = ChebyshevRate(Tmin, Tmax, Pmin, Pmax, coeffs) + + cdef _legacy_call(self, float T, float P): + cdef CxxChebyshevReaction2* r = self.cr2() cdef double logT = np.log(T) cdef double recipT = 1/T cdef double logP = np.log10(P) @@ -929,6 +1581,20 @@ cdef class ChebyshevReaction(Reaction): r.rate.update_C(&logP) return r.rate.updateRC(logT, recipT) + def __call__(self, float T, float P): + """ + .. deprecated:: 2.6 + To be deprecated with version 2.6, and removed thereafter. + Replaceable by call to `rate` property. + """ + if self.uses_legacy: + return self._legacy_call(T, P) + + warnings.warn( + self._deprecation_warning("__call__", "method"), DeprecationWarning) + return self.rate(T, P) + + cdef class BlowersMasel: """ A reaction rate coefficient which depends on temperature and enthalpy change @@ -996,18 +1662,20 @@ cdef class BlowersMasel: cdef double recipT = 1/T return self.rate.updateRC(logT, recipT, dH) + cdef wrapBlowersMasel(CxxBlowersMasel* rate, Reaction reaction): r = BlowersMasel(init=False) r.rate = rate r.reaction = reaction return r + cdef class BlowersMaselReaction(Reaction): """ A reaction which follows mass-action kinetics with Blowers Masel reaction rate. """ - reaction_type = "Blowers-Masel" + _reaction_type = "Blowers-Masel" property rate: """ Get/Set the `Arrhenius` rate coefficient for this reaction. """ @@ -1030,106 +1698,17 @@ cdef class BlowersMaselReaction(Reaction): cdef CxxBlowersMaselReaction* r = self.reaction r.allow_negative_pre_exponential_factor = allow -cdef class CustomReaction(Reaction): - """ - A reaction which follows mass-action kinetics with a custom reaction rate. - - An example for the definition of a `CustomReaction` object is given as:: - - rxn = CustomReaction(equation='H2 + O <=> H + OH', - rate=lambda T: 38.7 * T**2.7 * exp(-3150.15428/T), - kinetics=gas) - - Warning: this class is an experimental part of the Cantera API and - may be changed or removed without notice. - """ - reaction_type = "custom-rate-function" - - def __init__(self, equation=None, rate=None, Kinetics kinetics=None, - init=True, **kwargs): - - if init and equation and kinetics: - - yaml = '{{equation: {}, type: {}}}'.format(equation, self.reaction_type) - self._reaction = CxxNewReaction(AnyMapFromYamlString(stringify(yaml)), - deref(kinetics.kinetics)) - self.reaction = self._reaction.get() - self.rate = CustomRate(rate) - - property rate: - """ Get/Set the `CustomRate` object for this reaction. """ - def __get__(self): - return self._rate - def __set__(self, CustomRate rate): - self._rate = rate - cdef CxxCustomFunc1Reaction* r = self.reaction - r.setRate(self._rate._base) - - -cdef class TestReaction(Reaction): - """ - A reaction which follows mass-action kinetics with a modified Arrhenius - reaction rate. The class is a re-implementation of `ElementaryReaction` - and serves for testing purposes. - - An example for the definition of a `TestReaction` object is given as:: - - rxn = TestReaction(equation='H2 + O <=> H + OH', - rate={'A': 38.7, 'b': 2.7, 'Ea': 2.619184e+07}, - kinetics=gas) - - Warning: this class is an experimental part of the Cantera API and - may be changed or removed without notice. - """ - reaction_type = "elementary-new" - - def __init__(self, equation=None, rate=None, Kinetics kinetics=None, - init=True, **kwargs): - - if init and equation and kinetics: - - if isinstance(rate, dict): - coeffs = ['{}: {}'.format(k, v) for k, v in rate.items()] - elif isinstance(rate, ArrheniusRate) or rate is None: - coeffs = ['{}: 0.'.format(k) for k in ['A', 'b', 'Ea']] - else: - raise TypeError("Invalid rate definition") - - rate_def = '{{{}}}'.format(', '.join(coeffs)) - yaml = '{{equation: {}, rate-constant: {}, type: {}}}'.format( - equation, rate_def, self.reaction_type) - self._reaction = CxxNewReaction(AnyMapFromYamlString(stringify(yaml)), - deref(kinetics.kinetics)) - self.reaction = self._reaction.get() - - if isinstance(rate, ArrheniusRate): - self.rate = rate - - property rate: - """ Get/Set the `Arrhenius` rate coefficient for this reaction. """ - def __get__(self): - cdef CxxTestReaction* r = self.reaction - return ArrheniusRate.wrap(r.rate()) - def __set__(self, ArrheniusRate rate): - cdef CxxTestReaction* r = self.reaction - r.setRate(rate._base) - - property allow_negative_pre_exponential_factor: - """ - Get/Set whether the rate coefficient is allowed to have a negative - pre-exponential factor. - """ - def __get__(self): - cdef CxxElementaryReaction* r = self.reaction - return r.allow_negative_pre_exponential_factor - def __set__(self, allow): - cdef CxxElementaryReaction* r = self.reaction - r.allow_negative_pre_exponential_factor = allow - cdef class InterfaceReaction(ElementaryReaction): """ A reaction occurring on an `Interface` (i.e. a surface or an edge) """ - reaction_type = "interface" + _reaction_type = "interface" + _has_legacy = False + _hybrid = False + + property uses_legacy: + # legacy framework is implicitly used + def __get__(self): + return True property coverage_deps: """ @@ -1196,12 +1775,13 @@ cdef class InterfaceReaction(ElementaryReaction): cdef CxxInterfaceReaction* r = self.reaction r.sticking_species = stringify(species) + cdef class BlowersMaselInterfaceReaction(BlowersMaselReaction): """ A reaction occurring on an `Interface` (i.e. a surface or an edge) with the rate parameterization of `BlowersMasel`. """ - reaction_type = "surface-Blowers-Masel" + _reaction_type = "surface-Blowers-Masel" property coverage_deps: """ @@ -1267,3 +1847,41 @@ cdef class BlowersMaselInterfaceReaction(BlowersMaselReaction): def __set__(self, species): cdef CxxBlowersMaselInterfaceReaction* r = self.reaction r.sticking_species = stringify(species) + + +cdef class CustomReaction(Reaction): + """ + A reaction which follows mass-action kinetics with a custom reaction rate. + + An example for the definition of a `CustomReaction` object is given as:: + + rxn = CustomReaction( + equation="H2 + O <=> H + OH", + rate=lambda T: 38.7 * T**2.7 * exp(-3150.15428/T), + kinetics=gas) + + Warning: this class is an experimental part of the Cantera API and + may be changed or removed without notice. + """ + _reaction_type = "custom-rate-function" + + def __init__(self, equation=None, rate=None, Kinetics kinetics=None, + init=True, **kwargs): + + if init and equation and kinetics: + + spec = {"equation": equation, "type": self._reaction_type} + + self._reaction = CxxNewReaction(dict_to_anymap(spec), + deref(kinetics.kinetics)) + self.reaction = self._reaction.get() + self.rate = CustomRate(rate) + + property rate: + """ Get/Set the `CustomRate` object for this reaction. """ + def __get__(self): + return self._rate + def __set__(self, CustomRate rate): + self._rate = rate + cdef CxxCustomFunc1Reaction* r = self.reaction + r.setRate(self._rate._base) diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index a5e035ec8c1..db2becc7634 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -48,8 +48,8 @@ def test_multiplier(self): self.assertArrayNear(0.5 * rev_rates0, rev_rates2) def test_reaction_type(self): - self.assertEqual(self.phase.reaction_type_str(0), "three-body") - self.assertEqual(self.phase.reaction_type_str(2), "elementary") + self.assertIn(self.phase.reaction_type_str(0), ["three-body", "three-body-legacy"]) + self.assertIn(self.phase.reaction_type_str(2), ["elementary", "elementary-legacy"]) self.assertEqual(self.phase.reaction_type_str(21), "falloff") with self.assertRaisesRegex(ValueError, 'out of range'): @@ -404,7 +404,8 @@ def test_modify_thermo(self): self.assertAlmostEqual(w2[i] / w1[i], 1.0) def test_pdep_err(self): - err_msg = ("Invalid rate coefficient for reaction 'CH2CHOO <=> CH3O + CO'", + err_msg = ("CanteraError thrown by Plog::validate:", + "Invalid rate coefficient for reaction 'CH2CHOO <=> CH3O + CO'", "at P = 32019, T = 500.0", "at P = 32019, T = 1000.0", "at P = 1.0132e+05, T = 500.0", @@ -958,7 +959,7 @@ def test_input_data_from_file(self): def test_input_data_from_scratch(self): r = ct.ElementaryReaction({'O':1, 'H2':1}, {'H':1, 'OH':1}) - r.rate = ct.Arrhenius(3.87e1, 2.7, 2.6e7) + r.rate = ct.ArrheniusRate(3.87e1, 2.7, 2.6e7) data = r.input_data self.assertNear(data['rate-constant']['A'], 3.87e1) self.assertNear(data['rate-constant']['b'], 2.7) @@ -969,7 +970,7 @@ def test_input_data_from_scratch(self): def test_elementary(self): r = ct.ElementaryReaction({'O':1, 'H2':1}, {'H':1, 'OH':1}) - r.rate = ct.Arrhenius(3.87e1, 2.7, 6260*1000*4.184) + r.rate = ct.ArrheniusRate(3.87e1, 2.7, 6260*1000*4.184) gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=self.species, reactions=[r]) @@ -987,15 +988,15 @@ def test_arrhenius_rate(self): def test_negative_A(self): species = ct.Species.listFromFile('gri30.cti') r = ct.ElementaryReaction('NH:1, NO:1', 'N2O:1, H:1') - r.rate = ct.Arrhenius(-2.16e13, -0.23, 0) + r.rate = ct.ArrheniusRate(-2.16e13, -0.23, 0) - self.assertFalse(r.allow_negative_pre_exponential_factor) + self.assertFalse(r.rate.allow_negative_pre_exponential_factor) with self.assertRaisesRegex(ct.CanteraError, 'negative pre-exponential'): gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=species, reactions=[r]) - r.allow_negative_pre_exponential_factor = True + r.rate.allow_negative_pre_exponential_factor = True gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=species, reactions=[r]) @@ -1025,7 +1026,7 @@ def test_threebody(self): r = ct.ThreeBodyReaction() r.reactants = {'O':1, 'H':1} r.products = {'OH':1} - r.rate = ct.Arrhenius(5e11, -1.0, 0.0) + r.rate = ct.ArrheniusRate(5e11, -1.0, 0.0) r.efficiencies = {'AR':0.7, 'H2':2.0, 'H2O':6.0} gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', @@ -1061,7 +1062,7 @@ def test_plog(self): r = ct.PlogReaction() r.reactants = {'R1A':1, 'R1B':1} r.products = {'P1':1, 'H':1} - r.rates = [ + r.rate.rates = [ (0.01*ct.one_atm, ct.Arrhenius(1.2124e13, -0.5779, 10872.7*4184)), (1.0*ct.one_atm, ct.Arrhenius(4.9108e28, -4.8507, 24772.8*4184)), (10.0*ct.one_atm, ct.Arrhenius(1.2866e44, -9.0246, 39796.5*4184)), @@ -1084,7 +1085,7 @@ def test_plog_rate(self): gas1 = ct.Solution('pdep-test.yaml') gas1.TP = 800, 2*ct.one_atm for i in range(4): - self.assertNear(gas1.reaction(i)(gas1.T, gas1.P), + self.assertNear(gas1.reaction(i).rate(gas1.T, gas1.P), gas1.forward_rate_constants[i]) def test_chebyshev(self): @@ -1094,11 +1095,11 @@ def test_chebyshev(self): r = ct.ChebyshevReaction() r.reactants = 'R5:1, H:1' r.products = 'P5A:1, P5B:1' - r.set_parameters(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, - coeffs=[[ 5.28830e+00, -1.13970e+00, -1.20590e-01, 1.60340e-02], - [ 1.97640e+00, 1.00370e+00, 7.28650e-03, -3.04320e-02], - [ 3.17700e-01, 2.68890e-01, 9.48060e-02, -7.63850e-03], - [-3.12850e-02, -3.94120e-02, 4.43750e-02, 1.44580e-02]]) + r.rate = ct.ChebyshevRate(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, + data=[[ 5.28830e+00, -1.13970e+00, -1.20590e-01, 1.60340e-02], + [ 1.97640e+00, 1.00370e+00, 7.28650e-03, -3.04320e-02], + [ 3.17700e-01, 2.68890e-01, 9.48060e-02, -7.63850e-03], + [-3.12850e-02, -3.94120e-02, 4.43750e-02, 1.44580e-02]]) gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=species, reactions=[r]) @@ -1117,11 +1118,11 @@ def test_chebyshev_single_P(self): r = ct.ChebyshevReaction() r.reactants = 'R5:1, H:1' r.products = 'P5A:1, P5B:1' - r.set_parameters(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, - coeffs=[[ 5.28830e+00], - [ 1.97640e+00], - [ 3.17700e-01], - [-3.12850e-02]]) + r.rate = ct.ChebyshevRate(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, + data=[[ 5.28830e+00], + [ 1.97640e+00], + [ 3.17700e-01], + [-3.12850e-02]]) gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=species, reactions=[r]) @@ -1139,8 +1140,8 @@ def test_chebyshev_single_T(self): r = ct.ChebyshevReaction() r.reactants = 'R5:1, H:1' r.products = 'P5A:1, P5B:1' - r.set_parameters(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, - coeffs=[[ 5.28830e+00, -1.13970e+00, -1.20590e-01, 1.60340e-02]]) + r.rate = ct.ChebyshevRate(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, + data=[[ 5.28830e+00, -1.13970e+00, -1.20590e-01, 1.60340e-02]]) gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=species, reactions=[r]) @@ -1157,7 +1158,7 @@ def test_chebyshev_rate(self): gas1 = ct.Solution('pdep-test.yaml') gas1.TP = 800, 2*ct.one_atm for i in range(4,6): - self.assertNear(gas1.reaction(i)(gas1.T, gas1.P), + self.assertNear(gas1.reaction(i).rate(gas1.T, gas1.P), gas1.forward_rate_constants[i]) def test_chebyshev_bad_shape_cti(self): @@ -1349,7 +1350,7 @@ def test_modify_elementary(self): A2 = 1.5 * A1 b2 = b1 + 0.1 Ta2 = Ta1 * 1.2 - R.rate = ct.Arrhenius(A2, b2, Ta2 * ct.gas_constant) + R.rate = ct.ArrheniusRate(A2, b2, Ta2 * ct.gas_constant) gas.modify_reaction(2, R) self.assertNear(A2*T**b2*np.exp(-Ta2/T), gas.forward_rate_constants[2]) @@ -1364,7 +1365,7 @@ def test_modify_third_body(self): A2 = 1.7 * A1 b2 = b1 - 0.1 - R.rate = ct.Arrhenius(A2, b2, 0.0) + R.rate = ct.ArrheniusRate(A2, b2, 0.0) gas.modify_reaction(5, R) kf2 = gas.forward_rate_constants[5] self.assertNear((A2*T**b2) / (A1*T**b1), kf2/kf1) @@ -1393,13 +1394,13 @@ def test_modify_plog(self): r0 = gas.reaction(0) r1 = gas.reaction(1) - r0.rates = r1.rates + r0.rate = ct.PlogRate(r1.rate.rates) gas.modify_reaction(0, r0) kf = gas.forward_rate_constants self.assertNear(kf[0], kf[1]) # Removing the high-pressure rates should have no effect at low P... - r1.rates = r1.rates[:-4] + r1.rate = ct.PlogRate(rates=r1.rate.rates[:-4]) gas.modify_reaction(1, r1) self.assertNear(kf[1], gas.forward_rate_constants[1]) @@ -1414,7 +1415,8 @@ def test_modify_chebyshev(self): r1 = gas.reaction(4) r2 = gas.reaction(5) - r1.set_parameters(r2.Tmin, r2.Tmax, r2.Pmin, r2.Pmax, r2.coeffs) + r1.rate = ct.ChebyshevRate(r2.rate.Tmin, r2.rate.Tmax, + r2.rate.Pmin, r2.rate.Pmax, r2.rate.coeffs) # rates should be different before calling 'modify_reaction' kf = gas.forward_rate_constants diff --git a/interfaces/cython/cantera/test/test_reaction.py b/interfaces/cython/cantera/test/test_reaction.py index f2d66c6b36c..df21c31dfbf 100644 --- a/interfaces/cython/cantera/test/test_reaction.py +++ b/interfaces/cython/cantera/test/test_reaction.py @@ -2,10 +2,12 @@ from pathlib import Path import cantera as ct +import numpy as np from . import utilities class TestImplicitThirdBody(utilities.CanteraTest): + # tests for three-body reactions with specified collision partner @classmethod def setUpClass(cls): @@ -13,6 +15,7 @@ def setUpClass(cls): cls.gas = ct.Solution("gri30.yaml") def test_implicit_three_body(self): + # check equivalency of auto-detected and explicit specification yaml1 = """ equation: H + 2 O2 <=> HO2 + O2 rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} @@ -72,6 +75,7 @@ def test_duplicate(self): Path(fname).unlink() def test_short_serialization(self): + # check that serialized output is compact yaml = """ equation: H + O2 + H2O <=> HO2 + H2O rate-constant: {A: 1.126e+19, b: -0.76, Ea: 0.0} @@ -84,14 +88,16 @@ def test_short_serialization(self): self.assertNotIn("efficiencies", input_data) def test_non_integer_stoich(self): + # check that non-integer coefficients prevent automatic conversion yaml = """ - equation: H + 1.5 O2 <=> HO2 + O2 + equation: 2 H + 1.5 O2 <=> H2O + O2 rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} """ rxn = ct.Reaction.fromYaml(yaml, self.gas) self.assertEqual(rxn.reaction_type, "elementary") def test_not_three_body(self): + # check that insufficient reactants prevent automatic conversion yaml = """ equation: HCNO + H <=> H + HNCO # Reaction 270 rate-constant: {A: 2.1e+15, b: -0.69, Ea: 2850.0} @@ -100,6 +106,7 @@ def test_not_three_body(self): self.assertEqual(rxn.reaction_type, "elementary") def test_user_override(self): + # check that type specification prevents automatic conversion yaml = """ equation: H + 2 O2 <=> HO2 + O2 rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} @@ -109,78 +116,312 @@ def test_user_override(self): self.assertEqual(rxn.reaction_type, "elementary") -class TestElementary(utilities.CanteraTest): +class ReactionRateTests: + # test suite for reaction rate expressions - _cls = ct.ElementaryReaction - _equation = 'H2 + O <=> H + OH' - _rate = {'A': 38.7, 'b': 2.7, 'Ea': 2.619184e+07} - _rate_obj = ct.Arrhenius(38.7, 2.7, 2.619184e+07) - _index = 2 - _type = "elementary" + _cls = None # reaction rate object to be tested + _type = None # name of reaction rate + _uses_pressure = False # rate evaluation requires pressure + _index = None # index of reaction in "kineticsfromscratch.yaml" + _input = None # input parameters (dict corresponding to YAML) @classmethod def setUpClass(cls): utilities.CanteraTest.setUpClass() - cls.gas = ct.Solution('h2o2.yaml', transport_model=None) - cls.gas.X = 'H2:0.1, H2O:0.2, O2:0.7, O:1e-4, OH:1e-5, H:2e-5' + cls.gas = ct.Solution("kineticsfromscratch.yaml") + cls.gas.X = "H2:0.1, H2O:0.2, O2:0.7, O:1e-4, OH:1e-5, H:2e-5" cls.gas.TP = 900, 2*ct.one_atm - cls.species = ct.Species.listFromFile('h2o2.yaml') - def check_rxn(self, rxn): + def test_type(self): + # check reaction type + self.assertIn(self._type, "{}".format(self.rate)) + + def test_rate_T(self): + # check evaluation as a function of temperature only + if self._uses_pressure: + with self.assertRaisesRegex(ct.CanteraError, "reaction type requires pressure"): + self.assertNear(self.rate(self.gas.T), + self.gas.forward_rate_constants[self._index]) + else: + self.assertNear(self.rate(self.gas.T), + self.gas.forward_rate_constants[self._index]) + + def test_rate_TP(self): + # check evaluation as a function of temperature and pressure + self.assertNear(self.rate(self.gas.T, self.gas.P), + self.gas.forward_rate_constants[self._index]) + + def test_input(self): + # check instantiation based on input_data + rate = self._cls(input_data=self._input) + self.assertIn(self._type, "{}".format(rate)) + self.assertNear(rate(self.gas.T, self.gas.P), + self.rate(self.gas.T, self.gas.P)) + + def test_unconfigured(self): + # check behavior of unconfigured rate object + rate = self._cls(input_data={}) + self.assertTrue(np.isnan(rate(self.gas.T, self.gas.P))) + input_data = rate.input_data + self.assertIsInstance(input_data, dict) + self.assertEqual(input_data, {}) + + def test_roundtrip(self): + # check round-trip instantiation via input_data + if self._index is None: + return + input_data = self.rate.input_data + rate = self._cls(input_data=input_data) + self.assertNear(rate(self.gas.T, self.gas.P), + self.rate(self.gas.T, self.gas.P)) + + +class TestArrheniusRate(ReactionRateTests, utilities.CanteraTest): + # test Arrhenius rate expressions + + _cls = ct.ArrheniusRate + _type = "Arrhenius" + _uses_pressure = False + _index = 0 + _input = {"rate-constant": {"A": 38.7, "b": 2.7, "Ea": 26191840.0}} + + def setUp(self): + self.A = self.gas.reaction(self._index).rate.pre_exponential_factor + self.b = self.gas.reaction(self._index).rate.temperature_exponent + self.Ea = self.gas.reaction(self._index).rate.activation_energy + self.rate = ct.ArrheniusRate(self.A, self.b, self.Ea) + + def test_parameters(self): + # test reaction rate parameters + self.assertEqual(self.A, self.rate.pre_exponential_factor) + self.assertEqual(self.b, self.rate.temperature_exponent) + self.assertEqual(self.Ea, self.rate.activation_energy) + + def test_allow_negative_pre_exponential_factor(self): + # test reaction rate property + self.assertFalse(self.rate.allow_negative_pre_exponential_factor) + self.rate.allow_negative_pre_exponential_factor = True + self.assertTrue(self.rate.allow_negative_pre_exponential_factor) + + +class TestPlogRate(ReactionRateTests, utilities.CanteraTest): + # test Plog rate expressions + + _cls = ct.PlogRate + _type = "Plog" + _uses_pressure = True + _index = 3 + _input = {"rate-constants": [ + {"P": 1013.25, "A": 1.2124e+16, "b": -0.5779, "Ea": 45491376.8}, + {"P": 101325., "A": 4.9108e+31, "b": -4.8507, "Ea": 103649395.2}, + {"P": 1013250., "A": 1.2866e+47, "b": -9.0246, "Ea": 166508556.0}, + {"P": 10132500., "A": 5.9632e+56, "b": -11.529, "Ea": 220076726.4}]} + + def setUp(self): + self.rate = ct.PlogRate([(1013.25, ct.Arrhenius(1.2124e+16, -0.5779, 45491376.8)), + (101325., ct.Arrhenius(4.9108e+31, -4.8507, 103649395.2)), + (1013250., ct.Arrhenius(1.2866e+47, -9.0246, 166508556.0)), + (10132500., ct.Arrhenius(5.9632e+56, -11.529, 220076726.4))]) + + def test_get_rates(self): + # test getter for property rates + rates = self.rate.rates + self.assertIsInstance(rates, list) + + other = self._input["rate-constants"] + self.assertEqual(len(rates), len(other)) + for index, item in enumerate(rates): + P, rate = item + self.assertNear(P, other[index]["P"]) + self.assertNear(rate.pre_exponential_factor, other[index]["A"]) + self.assertNear(rate.temperature_exponent, other[index]["b"]) + self.assertNear(rate.activation_energy, other[index]["Ea"]) + + def test_set_rates(self): + # test setter for property rates + other = [ + {"P": 100., "A": 1.2124e+16, "b": -1., "Ea": 45491376.8}, + {"P": 10000., "A": 4.9108e+31, "b": -2., "Ea": 103649395.2}, + {"P": 1000000., "A": 1.2866e+47, "b": -3., "Ea": 166508556.0}] + rate = ct.PlogRate([(o["P"], ct.Arrhenius(o["A"], o["b"], o["Ea"])) + for o in other]) + rates = rate.rates + self.assertEqual(len(rates), len(other)) + + for index, item in enumerate(rates): + P, rate = item + self.assertNear(P, other[index]["P"]) + self.assertNear(rate.pre_exponential_factor, other[index]["A"]) + self.assertNear(rate.temperature_exponent, other[index]["b"]) + self.assertNear(rate.activation_energy, other[index]["Ea"]) + + def test_no_rates(self): + # test instantiation of empty rate + rate = ct.PlogRate() + self.assertIsInstance(rate.rates, list) + + +class TestChebyshevRate(ReactionRateTests, utilities.CanteraTest): + # test Chebyshev rate expressions + + _cls = ct.ChebyshevRate + _type = "Chebyshev" + _uses_pressure = True + _index = 4 + _input = {"data": [[8.2883, -1.1397, -0.12059, 0.016034], + [1.9764, 1.0037, 0.0072865, -0.030432], + [0.3177, 0.26889, 0.094806, -0.0076385]], + "pressure-range": [1000.0, 10000000.0], + "temperature-range": [290.0, 3000.0]} + + def setUp(self): + self.Tmin = self.gas.reaction(self._index).rate.Tmin + self.Tmax = self.gas.reaction(self._index).rate.Tmax + self.Pmin = self.gas.reaction(self._index).rate.Pmin + self.Pmax = self.gas.reaction(self._index).rate.Pmax + self.coeffs = self.gas.reaction(self._index).rate.coeffs + self.rate = ct.ChebyshevRate(self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.coeffs) + + def test_parameters(self): + # test getters for rate properties + self.assertEqual(self.Tmin, self.rate.Tmin) + self.assertEqual(self.Tmax, self.rate.Tmax) + self.assertEqual(self.Pmin, self.rate.Pmin) + self.assertEqual(self.Pmax, self.rate.Pmax) + self.assertTrue(np.all(self.coeffs == self.rate.coeffs)) + + +class ReactionTests: + # test suite for reaction expressions + + _cls = None # reaction object to be tested + _type = None # name of reaction rate + _legacy = False # object uses legacy framework + _equation = None # reaction equation string + _rate = None # parameters for reaction rate object constructor + _rate_obj = None # reaction rate object + _kwargs = {} # additional parameters required by contructor + _index = None # index of reaction in "kineticsfromscratch.yaml" + _yaml = None # YAML parameterization + _input = None # input parameters (dict corresponding to YAML) + _deprecated_getters = {} # test legacy getters (old framework) + _deprecated_setters = {} # test legacy setters (old framework) + _deprecated_callers = {} # test legacy callers (old framework) + + @classmethod + def setUpClass(cls): + utilities.CanteraTest.setUpClass() + cls.gas = ct.Solution("kineticsfromscratch.yaml", transport_model=None) + cls.gas.X = "H2:0.1, H2O:0.2, O2:0.7, O:1e-4, OH:1e-5, H:2e-5" + cls.gas.TP = 900, 2*ct.one_atm + cls.species = cls.gas.species() + + def check_rxn(self, rxn, check_legacy=True): + # helper function that checks reaction configuration ix = self._index - self.assertEqual(rxn.reaction_type, self._type) - self.assertNear(rxn.rate(self.gas.T), self.gas.forward_rate_constants[ix]) self.assertEqual(rxn.reactants, self.gas.reaction(ix).reactants) self.assertEqual(rxn.products, self.gas.reaction(ix).products) + if check_legacy: + self.assertEqual(rxn.reaction_type, self._type) + self.assertEqual(rxn.uses_legacy, self._type.endswith("-legacy")) + self.assertEqual(rxn.uses_legacy, self._legacy) - gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + gas2 = ct.Solution(thermo="IdealGas", kinetics="GasKinetics", species=self.species, reactions=[rxn]) gas2.TPX = self.gas.TPX - self.check_sol(gas2) + self.check_solution(gas2, check_legacy) - def check_sol(self, gas2): + def check_solution(self, gas2, check_legacy=True): + # helper function that checks evaluation of reaction rates ix = self._index - self.assertEqual(gas2.reaction_type_str(0), self._type) + if check_legacy: + self.assertEqual(gas2.reaction_type_str(0), self._type) self.assertNear(gas2.forward_rate_constants[0], self.gas.forward_rate_constants[ix]) self.assertNear(gas2.net_rates_of_progress[0], self.gas.net_rates_of_progress[ix]) def test_rate(self): - self.assertNear(self._rate_obj(self.gas.T), self.gas.forward_rate_constants[self._index]) + # check consistency of reaction rate and forward rate constant + if self._rate_obj is None: + return + if self._legacy: + self.assertNear(self._rate_obj(self.gas.T), self.gas.forward_rate_constants[self._index]) + else: + self.assertNear(self._rate_obj(self.gas.T, self.gas.P), + self.gas.forward_rate_constants[self._index]) def test_from_parts(self): + # check instantiation from parts (reactants, products, rate expression) + if not self._rate_obj: + return orig = self.gas.reaction(self._index) - rxn = self._cls(orig.reactants, orig.products) + rxn = self._cls(orig.reactants, orig.products, legacy=self._legacy) rxn.rate = self._rate_obj self.check_rxn(rxn) - def test_from_dict(self): - rxn = self._cls(equation=self._equation, rate=self._rate, kinetics=self.gas) + def test_from_dict1(self): + # check instantiation from keywords / rate defined by dictionary + rxn = self._cls(equation=self._equation, rate=self._rate, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) self.check_rxn(rxn) + def test_from_yaml(self): + # check instantiation from yaml string + if self._yaml is None: + return + rxn = ct.Reaction.fromYaml(self._yaml, kinetics=self.gas) + self.check_rxn(rxn) + + def test_from_dict2(self): + # check instantiation from a yaml dictionary + if self._yaml is None: + return + rxn1 = ct.Reaction.fromYaml(self._yaml, kinetics=self.gas) + input_data = rxn1.input_data + rxn2 = ct.Reaction.from_dict(input_data, kinetics=self.gas) + # cannot compare types as input_data does not recreate legacy objects + self.check_rxn(rxn2, check_legacy=False) + def test_from_rate(self): - rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas) + # check instantiation from keywords / rate provided as object + if self._rate_obj is None: + return + rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) self.check_rxn(rxn) def test_add_rxn(self): - gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + # check adding new reaction to solution + if self._rate_obj is None: + return + gas2 = ct.Solution(thermo="IdealGas", kinetics="GasKinetics", species=self.species, reactions=[]) gas2.TPX = self.gas.TPX - rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas) + rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) gas2.add_reaction(rxn) - self.check_sol(gas2) + self.check_solution(gas2) - def test_wrong_rate(self): + def test_raises_invalid_rate(self): + # check exception for instantiation from keywords / invalid rate with self.assertRaises(TypeError): - rxn = self._cls(equation=self._equation, rate=[], kinetics=self.gas) + rxn = self._cls(equation=self._equation, rate=(), kinetics=self.gas, + legacy=self._legacy, **self._kwargs) def test_no_rate(self): - rxn = self._cls(equation=self._equation, kinetics=self.gas) - self.assertNear(rxn.rate(self.gas.T), 0.) - - gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + # check behavior for instantiation from keywords / no rate + if self._rate_obj is None: + return + rxn = self._cls(equation=self._equation, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) + if self._legacy: + self.assertNear(rxn.rate(self.gas.T), 0.) + else: + self.assertTrue(np.isnan(rxn.rate(self.gas.T, self.gas.P))) + + gas2 = ct.Solution(thermo="IdealGas", kinetics="GasKinetics", species=self.species, reactions=[rxn]) gas2.TPX = self.gas.TPX @@ -188,58 +429,378 @@ def test_no_rate(self): self.assertNear(gas2.net_rates_of_progress[0], 0.) def test_replace_rate(self): - rxn = self._cls(equation=self._equation, kinetics=self.gas) + # check replacing reaction rate expression + if self._rate_obj is None: + return + rxn = self._cls(equation=self._equation, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) + rxn.rate = self._rate_obj + self.check_rxn(rxn) + + def test_roundtrip(self): + # check round-trip instantiation via input_data + if self._legacy: + return + rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) + input_data = rxn.rate.input_data + rate_obj = rxn.rate.__class__(input_data=input_data) + rxn2 = self._cls(equation=self._equation, rate=rate_obj, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) + self.check_rxn(rxn2) + + def check_equal(self, one, two): + # helper function for deprecation tests + self.assertEqual(type(one), type(two)) + if isinstance(one, (list, tuple, np.ndarray)): + self.assertArrayNear(one, two) + else: + self.assertEqual(one, two) + + def test_deprecated_getters(self): + # check property getters deprecated in new framework + if self._yaml is None: + return + + rxn = ct.Reaction.fromYaml(self._yaml, kinetics=self.gas) + for attr, default in self._deprecated_getters.items(): + if self._legacy: + self.check_equal(getattr(rxn, attr), default) + else: + with self.assertWarnsRegex(DeprecationWarning, "property is moved"): + self.check_equal(getattr(rxn, attr), default) + + def test_deprecated_setters(self): + # check property setters deprecated in new framework + if self._yaml is None: + return + + rxn = ct.Reaction.fromYaml(self._yaml, kinetics=self.gas) + for attr, new in self._deprecated_setters.items(): + if self._legacy: + setattr(rxn, attr, new) + self.check_equal(getattr(rxn, attr), new) + else: + with self.assertWarnsRegex(DeprecationWarning, "property is moved"): + setattr(rxn, attr, new) + with self.assertWarnsRegex(DeprecationWarning, "property is moved"): + self.check_equal(getattr(rxn, attr), new) + + def test_deprecated_callers(self): + # check methods deprecated in new framework + if self._yaml is None: + return + + rxn = ct.Reaction.fromYaml(self._yaml, kinetics=self.gas) + for state, value in self._deprecated_callers.items(): + T, P = state + if self._legacy: + self.check_equal(rxn(T, P), value) + else: + with self.assertWarnsRegex(DeprecationWarning, "method is moved"): + self.check_equal(rxn(T, P), value) + + +class TestElementary(ReactionTests, utilities.CanteraTest): + # test updated version of elementary reaction + + _cls = ct.ElementaryReaction + _type = "elementary" + _equation = "H2 + O <=> H + OH" + _rate = {"A": 38.7, "b": 2.7, "Ea": 2.619184e+07} + _index = 0 + _yaml = """ + equation: O + H2 <=> H + OH + rate-constant: {A: 38.7, b: 2.7, Ea: 6260.0 cal/mol} + """ + _deprecated_getters = {"allow_negative_pre_exponential_factor": False} + _deprecated_setters = {"allow_negative_pre_exponential_factor": True} + + @classmethod + def setUpClass(cls): + ReactionTests.setUpClass() + if cls._legacy: + args = list(cls._rate.values()) + cls._rate_obj = ct.Arrhenius(*args) + else: + cls._rate_obj = ct.ArrheniusRate(**cls._rate) + + def test_arrhenius(self): + # test assigning Arrhenius rate + rate = ct.Arrhenius(self._rate["A"], self._rate["b"], self._rate["Ea"]) + rxn = self._cls(equation=self._equation, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) + if self._legacy: + rxn.rate = rate + else: + with self.assertWarnsRegex(DeprecationWarning, "'Arrhenius' object is deprecated"): + rxn.rate = rate + self.check_rxn(rxn) + + +class TestElementary2(TestElementary): + # test legacy version of elementary reaction + + _cls = ct.ElementaryReaction + _type = "elementary-legacy" + _legacy = True + _yaml = """ + equation: O + H2 <=> H + OH + type: elementary-legacy + rate-constant: {A: 38.7, b: 2.7, Ea: 6260.0 cal/mol} + """ + + +class TestThreeBody(TestElementary): + # test updated version of three-body reaction + + _cls = ct.ThreeBodyReaction + _type = "three-body" + _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}} + _index = 1 + _yaml = """ + equation: 2 O + M <=> O2 + M + type: three-body + 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 test_from_parts(self): + # overload default reaction creation from parts + orig = self.gas.reaction(self._index) + rxn = self._cls(orig.reactants, orig.products, legacy=self._legacy) rxn.rate = self._rate_obj + rxn.efficiencies = self._kwargs["efficiencies"] self.check_rxn(rxn) + def test_rate(self): + # rate constant contains third-body concentration + pass + + def test_efficiencies(self): + # check efficiencies + rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas, + legacy=self._legacy, **self._kwargs) + + self.assertEqual(rxn.efficiencies, self._kwargs["efficiencies"]) + + +class TestThreeBody2(TestThreeBody): + # test legacy version of three-body reaction + + _cls = ct.ThreeBodyReaction + _type = "three-body-legacy" + _legacy = True + _yaml = """ + equation: 2 O + M <=> O2 + M + type: three-body-legacy + rate-constant: {A: 1.2e+11, b: -1.0, Ea: 0.0 cal/mol} + efficiencies: {H2: 2.4, H2O: 15.4, AR: 0.83} + """ + + +class TestImplicitThreeBody(TestThreeBody): + # test three-body reactions with explicit collision parther + + _cls = ct.ThreeBodyReaction + _type = "three-body" + _equation = "H + 2 O2 <=> HO2 + O2" + _rate = {"A": 2.08e+19, "b": -1.24, "Ea": 0.0} + _index = 5 + _yaml = """ + equation: H + 2 O2 <=> HO2 + O2 + rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} + """ + + def test_efficiencies(self): + # overload of default tester + rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas, + legacy=self._legacy) + self.assertEqual(rxn.efficiencies, {"O2": 1.}) + self.assertEqual(rxn.default_efficiency, 0.) + + def test_from_parts(self): + # overload of default tester + orig = self.gas.reaction(self._index) + rxn = self._cls(orig.reactants, orig.products) + rxn.rate = self._rate_obj + rxn.efficiencies = {"O2": 1.} + rxn.default_efficiency = 0 + self.check_rxn(rxn) + + +class TestPlog(ReactionTests, utilities.CanteraTest): + # test updated version of Plog reaction + + _cls = ct.PlogReaction + _type = "pressure-dependent-Arrhenius" + _legacy = False + _equation = "H2 + O2 <=> 2 OH" + _rate = [(1013.25, ct.Arrhenius(1.2124e+16, -0.5779, 45491376.8)), + (101325., ct.Arrhenius(4.9108e+31, -4.8507, 103649395.2)), + (1013250., ct.Arrhenius(1.2866e+47, -9.0246, 166508556.0)), + (10132500., ct.Arrhenius(5.9632e+56, -11.529, 220076726.4))] + _index = 3 + _yaml = """ + equation: H2 + O2 <=> 2 OH + type: pressure-dependent-Arrhenius + rate-constants: + - {P: 0.01 atm, A: 1.2124e+16, b: -0.5779, Ea: 1.08727e+04 cal/mol} + - {P: 1.0 atm, A: 4.9108e+31, b: -4.8507, Ea: 2.47728e+04 cal/mol} + - {P: 10.0 atm, A: 1.2866e+47, b: -9.0246, Ea: 3.97965e+04 cal/mol} + - {P: 100.0 atm, A: 5.9632e+56, b: -11.529, Ea: 5.25996e+04 cal/mol} + """ + _deprecated_callers = {(1000., ct.one_atm): 530968934612.9017} -class TestCustom(TestElementary): + @classmethod + def setUpClass(cls): + ReactionTests.setUpClass() + if not cls._legacy: + cls._rate_obj = ct.PlogRate(cls._rate) + + def check_rates(self, rates, other): + # helper function used by deprecation tests + self.assertEqual(len(rates), len(other)) + for index, item in enumerate(rates): + P, rate = item + self.assertNear(P, other[index][0]) + self.assertNear(rate.pre_exponential_factor, other[index][1].pre_exponential_factor) + self.assertNear(rate.temperature_exponent, other[index][1].temperature_exponent) + self.assertNear(rate.activation_energy, other[index][1].activation_energy) + + def test_deprecated_getters(self): + # overload default tester for deprecated property getters + rxn = ct.Reaction.fromYaml(self._yaml, kinetics=self.gas) + if self._legacy: + self.check_rates(rxn.rates, self._rate) + else: + with self.assertWarnsRegex(DeprecationWarning, "property is moved"): + self.check_rates(rxn.rates, self._rate) + + def test_deprecated_setters(self): + # overload default tester for deprecated property setters + rate = ct.PlogRate(self._rate) + rates = rate.rates + + rxn = ct.Reaction.fromYaml(self._yaml, kinetics=self.gas) + if self._legacy: + rxn.rates = rates + self.check_rates(rxn.rates, self._rate) + else: + with self.assertWarnsRegex(DeprecationWarning, "Setter is replaceable"): + rxn.rates = rates + with self.assertWarnsRegex(DeprecationWarning, "property is moved"): + self.check_rates(rxn.rates, self._rate) + + +class TestPlog2(TestPlog): + # test legacy version of Plog reaction + + _cls = ct.PlogReaction + _type = "pressure-dependent-Arrhenius-legacy" + _legacy = True + _rate_obj = None + _yaml = """ + equation: H2 + O2 <=> 2 OH + type: pressure-dependent-Arrhenius-legacy + rate-constants: + - {P: 0.01 atm, A: 1.2124e+16, b: -0.5779, Ea: 1.08727e+04 cal/mol} + - {P: 1.0 atm, A: 4.9108e+31, b: -4.8507, Ea: 2.47728e+04 cal/mol} + - {P: 10.0 atm, A: 1.2866e+47, b: -9.0246, Ea: 3.97965e+04 cal/mol} + - {P: 100.0 atm, A: 5.9632e+56, b: -11.529, Ea: 5.25996e+04 cal/mol} + """ + + +class TestChebyshev(ReactionTests, utilities.CanteraTest): + # test updated version of Chebyshev reaction + + _cls = ct.ChebyshevReaction + _type = "Chebyshev" + _equation = "HO2 <=> OH + O" + _rate = {"Tmin": 290., "Tmax": 3000., "Pmin": 1000., "Pmax": 10000000.0, + "data": [[ 8.2883e+00, -1.1397e+00, -1.2059e-01, 1.6034e-02], + [ 1.9764e+00, 1.0037e+00, 7.2865e-03, -3.0432e-02], + [ 3.1770e-01, 2.6889e-01, 9.4806e-02, -7.6385e-03]]} + _index = 4 + _yaml = """ + equation: HO2 <=> OH + O + type: Chebyshev + temperature-range: [290.0, 3000.0] + pressure-range: [9.869232667160128e-03 atm, 98.69232667160128 atm] + 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] + """ + _deprecated_getters = {"nPressure": 4, "nTemperature": 3} + _deprecated_callers = {(1000., ct.one_atm): 2858762454.1119065} + + @classmethod + def setUpClass(cls): + ReactionTests.setUpClass() + if not cls._legacy: + cls._rate_obj = ct.ChebyshevRate(**cls._rate) + cls._deprecated_getters.update({"coeffs": np.array(cls._rate["data"])}) + cls._deprecated_getters.update( + {k: v for k, v in cls._rate.items() if k != "data"}) + + +class TestChebyshev2(TestChebyshev): + # test legacy version of Chebyshev reaction + + _cls = ct.ChebyshevReaction + _type = "Chebyshev-legacy" + _legacy = True + _rate_obj = None + _yaml = """ + equation: HO2 <=> OH + O + type: Chebyshev-legacy + temperature-range: [290.0, 3000.0] + pressure-range: [9.869232667160128e-03 atm, 98.69232667160128 atm] + 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] + """ + + +class TestCustom(ReactionTests, utilities.CanteraTest): + # test Custom reaction # probe O + H2 <=> H + OH _cls = ct.CustomReaction - _equation = 'H2 + O <=> H + OH' - _rate_obj = ct.CustomRate(lambda T: 38.7 * T**2.7 * exp(-3150.15428/T)) - _index = 2 _type = "custom-rate-function" + _legacy = False + _equation = "H2 + O <=> H + OH" + _rate_obj = ct.CustomRate(lambda T: 38.7 * T**2.7 * exp(-3150.15428/T)) + _index = 0 + _yaml = None def setUp(self): - # need to overwrite rate to ensure correct type ('method' is not compatible with Func1) + # need to overwrite rate to ensure correct type ("method" is not compatible with Func1) self._rate = lambda T: 38.7 * T**2.7 * exp(-3150.15428/T) - def test_no_rate(self): - rxn = self._cls(equation=self._equation, kinetics=self.gas) - with self.assertRaisesRegex(ct.CanteraError, "Custom rate function is not initialized."): - rxn.rate(self.gas.T) + def test_roundtrip(self): + # overload default tester for round trip + pass - gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', - species=self.species, reactions=[rxn]) - gas2.TPX = self.gas.TPX - - with self.assertRaisesRegex(ct.CanteraError, "Custom rate function is not initialized."): - gas2.forward_rate_constants - - def test_from_func(self): + def test_from_func1(self): + # check instantiation from keywords / rate provided as func1 f = ct.Func1(self._rate) rxn = ct.CustomReaction(equation=self._equation, rate=f, kinetics=self.gas) self.check_rxn(rxn) def test_rate_func(self): + # check result of rate expression f = ct.Func1(self._rate) rate = ct.CustomRate(f) self.assertNear(rate(self.gas.T), self.gas.forward_rate_constants[self._index]) - def test_custom(self): + def test_custom_lambda(self): + # check instantiation from keywords / rate provided as lambda function rxn = ct.CustomReaction(equation=self._equation, rate=lambda T: 38.7 * T**2.7 * exp(-3150.15428/T), kinetics=self.gas) self.check_rxn(rxn) - - -class TestElementaryNew(TestElementary): - - _cls = ct.TestReaction - _equation = 'H2 + O <=> H + OH' - _rate = {'A': 38.7, 'b': 2.7, 'Ea': 2.619184e+07} - _rate_obj = ct.ArrheniusRate(38.7, 2.7, 2.619184e+07) - _index = 2 - _type = "elementary-new" diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index 270545709b9..36914c73f84 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -103,6 +103,7 @@ bool BulkKinetics::addReaction(shared_ptr r) { bool added = Kinetics::addReaction(r); if (!added) { + // undeclared species, etc. return false; } double dn = 0.0; @@ -121,9 +122,8 @@ bool BulkKinetics::addReaction(shared_ptr r) m_irrev.push_back(nReactions()-1); } - if (std::dynamic_pointer_cast(r) != nullptr) { - shared_ptr rate; - rate = std::dynamic_pointer_cast(r)->rate(); + if (!(r->usesLegacy())) { + shared_ptr rate = r->rate(); // If neccessary, add new MultiBulkRates evaluator if (m_bulk_types.find(rate->type()) == m_bulk_types.end()) { m_bulk_types[rate->type()] = m_bulk_rates.size(); @@ -131,21 +131,57 @@ bool BulkKinetics::addReaction(shared_ptr r) if (rate->type() == "ArrheniusRate") { m_bulk_rates.push_back(std::unique_ptr( new MultiBulkRates)); + } else if (rate->type() == "PlogRate") { + m_bulk_rates.push_back(std::unique_ptr( + new MultiBulkRates)); + } else if (rate->type() == "ChebyshevRate") { + m_bulk_rates.push_back(std::unique_ptr( + new MultiBulkRates)); } else if (rate->type() == "custom-function") { m_bulk_rates.push_back(std::unique_ptr( new MultiBulkRates)); + } else { + throw CanteraError("BulkKinetics::addReaction", "Adding " + "reaction type '" + rate->type() + "' is not implemented"); } + m_bulk_rates.back()->resizeSpecies(m_kk); } // Add reaction rate to evaluator size_t index = m_bulk_types[rate->type()]; m_bulk_rates[index]->add(nReactions() - 1, *rate); + + // Add reaction to third-body evaluator + if (r->thirdBody() != nullptr) { + addThirdBody(r); + } } + m_concm.push_back(0.0); + return true; } -void BulkKinetics::addElementaryReaction(ElementaryReaction& r) +void BulkKinetics::addThirdBody(shared_ptr r) +{ + std::map efficiencies; + for (const auto& eff : r->thirdBody()->efficiencies) { + size_t k = kineticsSpeciesIndex(eff.first); + if (k != npos) { + efficiencies[k] = eff.second; + } else if (!m_skipUndeclaredThirdBodies) { + throw CanteraError("BulkKinetics::addThirdBody", "Found " + "third-body efficiency for undefined species '" + eff.first + + "' while adding reaction '" + r->equation() + "'"); + } + } + m_multi_concm.install(nReactions() - 1, efficiencies, + r->thirdBody()->default_efficiency); + concm_multi_values.resize(m_multi_concm.workSize()); + m_multi_indices.push_back(nReactions() - 1); +} + +void BulkKinetics::addElementaryReaction(ElementaryReaction2& r) { m_rates.install(nReactions()-1, r.rate); } @@ -155,11 +191,10 @@ void BulkKinetics::modifyReaction(size_t i, shared_ptr rNew) // operations common to all reaction types Kinetics::modifyReaction(i, rNew); - if (std::dynamic_pointer_cast(rNew) != nullptr) { - shared_ptr rate; - rate = std::dynamic_pointer_cast(rNew)->rate(); + if (!(rNew->usesLegacy())) { + shared_ptr rate = rNew->rate(); // Ensure that MultiBulkRates evaluator is available - if (m_bulk_types.find(rate->type()) != m_bulk_types.end()) { + if (m_bulk_types.find(rate->type()) == m_bulk_types.end()) { throw CanteraError("BulkKinetics::modifyReaction", "Evaluator not available for type '{}'.", rate->type()); } @@ -170,7 +205,7 @@ void BulkKinetics::modifyReaction(size_t i, shared_ptr rNew) } } -void BulkKinetics::modifyElementaryReaction(size_t i, ElementaryReaction& rNew) +void BulkKinetics::modifyElementaryReaction(size_t i, ElementaryReaction2& rNew) { m_rates.replace(i, rNew.rate); } @@ -181,6 +216,9 @@ void BulkKinetics::resizeSpecies() m_act_conc.resize(m_kk); m_phys_conc.resize(m_kk); m_grt.resize(m_kk); + for (auto& rates : m_bulk_rates) { + rates->resizeSpecies(m_kk); + } } void BulkKinetics::setMultiplier(size_t i, double f) { diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index b89505d2b5b..8812c093a95 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -53,8 +53,8 @@ void GasKinetics::update_rates_T() // loop over MultiBulkRates evaluators for (auto& rates : m_bulk_rates) { - rates->update(thermo()); - rates->getRateConstants(thermo(), m_rfn.data()); + rates->update(thermo(), m_concm.data()); + rates->getRateConstants(thermo(), m_rfn.data(), m_concm.data()); } if (m_plog_rates.nReactions()) { @@ -87,6 +87,15 @@ void GasKinetics::update_rates_C() m_falloff_concm.update(m_phys_conc, ctot, concm_falloff_values.data()); } + // Third-body objects interacting with MultiRate evaluator + if (!concm_multi_values.empty()) { + // using pre-existing third-body handlers requires copying; + m_multi_concm.update(m_phys_conc, ctot, concm_multi_values.data()); + for (size_t i = 0; i < m_multi_indices.size(); i++) { + m_concm[m_multi_indices[i]] = concm_multi_values[i]; + } + } + // P-log reactions if (m_plog_rates.nReactions()) { double logP = log(thermo().pressure()); @@ -184,6 +193,11 @@ void GasKinetics::updateROP() processFalloffReactions(); } + // reactions involving third body + for (auto& index : m_multi_indices) { + m_ropf[index] *= m_concm[index]; + } + for (size_t i = 0; i < nReactions(); i++) { // Scale the forward rate coefficient by the perturbation factor m_ropf[i] *= m_perturb[i]; @@ -230,6 +244,11 @@ void GasKinetics::getFwdRateConstants(doublereal* kfwd) processFalloffReactions(); } + // reactions involving third body + for (auto& index : m_multi_indices) { + m_ropf[index] *= m_concm[index]; + } + for (size_t i = 0; i < nReactions(); i++) { // multiply by perturbation factor kfwd[i] = m_ropf[i] * m_perturb[i]; @@ -242,23 +261,23 @@ bool GasKinetics::addReaction(shared_ptr r) bool added = BulkKinetics::addReaction(r); if (!added) { return false; - } else if (std::dynamic_pointer_cast(r) != nullptr) { + } else if (!(r->usesLegacy())) { // Rate object already added in BulkKinetics::addReaction return true; } - if (r->type() == "elementary") { - addElementaryReaction(dynamic_cast(*r)); - } else if (r->type() == "three-body") { - addThreeBodyReaction(dynamic_cast(*r)); + if (r->type() == "elementary-legacy") { + addElementaryReaction(dynamic_cast(*r)); + } else if (r->type() == "three-body-legacy") { + addThreeBodyReaction(dynamic_cast(*r)); } else if (r->type() == "falloff") { addFalloffReaction(dynamic_cast(*r)); } else if (r->type() == "chemically-activated") { addFalloffReaction(dynamic_cast(*r)); - } else if (r->type() == "pressure-dependent-Arrhenius") { - addPlogReaction(dynamic_cast(*r)); - } else if (r->type() == "Chebyshev") { - addChebyshevReaction(dynamic_cast(*r)); + } else if (r->type() == "pressure-dependent-Arrhenius-legacy") { + addPlogReaction(dynamic_cast(*r)); + } else if (r->type() == "Chebyshev-legacy") { + addChebyshevReaction(dynamic_cast(*r)); } else if (r->type() == "Blowers-Masel") { addBlowersMaselReaction(dynamic_cast(*r)); } else { @@ -299,7 +318,7 @@ void GasKinetics::addFalloffReaction(FalloffReaction& r) falloff_work.resize(m_falloffn.workSize()); } -void GasKinetics::addThreeBodyReaction(ThreeBodyReaction& r) +void GasKinetics::addThreeBodyReaction(ThreeBodyReaction2& r) { m_rates.install(nReactions()-1, r.rate); map efficiencies; @@ -314,12 +333,12 @@ void GasKinetics::addThreeBodyReaction(ThreeBodyReaction& r) concm_3b_values.resize(m_3b_concm.workSize()); } -void GasKinetics::addPlogReaction(PlogReaction& r) +void GasKinetics::addPlogReaction(PlogReaction2& r) { m_plog_rates.install(nReactions()-1, r.rate); } -void GasKinetics::addChebyshevReaction(ChebyshevReaction& r) +void GasKinetics::addChebyshevReaction(ChebyshevReaction2& r) { m_cheb_rates.install(nReactions()-1, r.rate); } @@ -334,23 +353,23 @@ void GasKinetics::modifyReaction(size_t i, shared_ptr rNew) // operations common to all bulk reaction types BulkKinetics::modifyReaction(i, rNew); - if (std::dynamic_pointer_cast(rNew) != nullptr) { + if (!(rNew->usesLegacy())) { // Rate object already modified in BulkKinetics::modifyReaction return; } - if (rNew->type() == "elementary") { - modifyElementaryReaction(i, dynamic_cast(*rNew)); - } else if (rNew->type() == "three-body") { - modifyThreeBodyReaction(i, dynamic_cast(*rNew)); + if (rNew->type() == "elementary-legacy") { + modifyElementaryReaction(i, dynamic_cast(*rNew)); + } else if (rNew->type() == "three-body-legacy") { + modifyThreeBodyReaction(i, dynamic_cast(*rNew)); } else if (rNew->type() == "falloff") { modifyFalloffReaction(i, dynamic_cast(*rNew)); } else if (rNew->type() == "chemically-activated") { modifyFalloffReaction(i, dynamic_cast(*rNew)); - } else if (rNew->type() == "pressure-dependent-Arrhenius") { - modifyPlogReaction(i, dynamic_cast(*rNew)); - } else if (rNew->type() == "Chebyshev") { - modifyChebyshevReaction(i, dynamic_cast(*rNew)); + } else if (rNew->type() == "pressure-dependent-Arrhenius-legacy") { + modifyPlogReaction(i, dynamic_cast(*rNew)); + } else if (rNew->type() == "Chebyshev-legacy") { + modifyChebyshevReaction(i, dynamic_cast(*rNew)); } else if (rNew->type() == "Blowers-Masel") { modifyBlowersMaselReaction(i, dynamic_cast(*rNew)); } else { @@ -364,7 +383,7 @@ void GasKinetics::modifyReaction(size_t i, shared_ptr rNew) m_pres += 0.1234; } -void GasKinetics::modifyThreeBodyReaction(size_t i, ThreeBodyReaction& r) +void GasKinetics::modifyThreeBodyReaction(size_t i, ThreeBodyReaction2& r) { m_rates.replace(i, r.rate); } @@ -377,12 +396,12 @@ void GasKinetics::modifyFalloffReaction(size_t i, FalloffReaction& r) m_falloffn.replace(iFall, r.falloff); } -void GasKinetics::modifyPlogReaction(size_t i, PlogReaction& r) +void GasKinetics::modifyPlogReaction(size_t i, PlogReaction2& r) { m_plog_rates.replace(i, r.rate); } -void GasKinetics::modifyChebyshevReaction(size_t i, ChebyshevReaction& r) +void GasKinetics::modifyChebyshevReaction(size_t i, ChebyshevReaction2& r) { m_cheb_rates.replace(i, r.rate); } diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 4874d7d70ec..92b13370e2f 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -144,8 +144,24 @@ std::pair Kinetics::checkDuplicates(bool throw_err) const continue; // No overlap in third body efficiencies } } else if (R.type() == "three-body") { - ThirdBody& tb1 = dynamic_cast(R).third_body; - ThirdBody& tb2 = dynamic_cast(other).third_body; + ThirdBody& tb1 = *(dynamic_cast(R).thirdBody()); + ThirdBody& 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-legacy") { + ThirdBody& tb1 = dynamic_cast(R).third_body; + ThirdBody& tb2 = dynamic_cast(other).third_body; bool thirdBodyOk = true; for (size_t k = 0; k < nTotalSpecies(); k++) { string s = kineticsSpeciesName(k); @@ -338,7 +354,7 @@ double Kinetics::productStoichCoeff(size_t kSpec, size_t irxn) const } int Kinetics::reactionType(size_t i) const { - warn_deprecated("Kinetics::reactionType()", + warn_deprecated("Kinetics::reactionType", "To be changed after Cantera 2.6. " "Return string instead of magic number; use " "Kinetics::reactionTypeStr during transition."); diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 7e783631e2c..21eee166913 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -111,6 +111,11 @@ void Reaction::validate() throw InputFileError("Reaction::validate", input, "Reaction orders may only be given for irreversible reactions"); } + + // Call validation of reaction rate evaluator + if (!usesLegacy()) { + m_rate->validate(equation()); + } } AnyMap Reaction::parameters(bool withInput) const @@ -153,6 +158,38 @@ void Reaction::getParameters(AnyMap& reactionNode) const } } +void Reaction::setParameters(const AnyMap& node, const Kinetics& kin) +{ + if (node.empty()) { + // empty node: used by newReaction() factory loader + return; + } + + parseReactionEquation(*this, node["equation"], kin); + // Non-stoichiometric reaction orders + if (node.hasKey("orders")) { + for (const auto& order : node["orders"].asMap()) { + orders[order.first] = order.second; + if (kin.kineticsSpeciesIndex(order.first) == npos) { + setValid(false); + } + } + } + + // remove optional third body notation (used by ChebyshevReaction) + reactants.erase("(+M)"); + products.erase("(+M)"); + + // Flags + id = node.getString("id", ""); + duplicate = node.getBool("duplicate", false); + allow_negative_orders = node.getBool("negative-orders", false); + allow_nonreactant_orders = node.getBool("nonreactant-orders", false); + + calculateRateCoeffUnits(kin); + input = node; +} + std::string Reaction::reactantString() const { std::ostringstream result; @@ -236,6 +273,12 @@ std::pair, bool> Reaction::undeclaredThirdBodies( const Kinetics& kin) const { std::vector undeclared; + if (m_third_body) { + updateUndeclared(undeclared, m_third_body->efficiencies, kin); + bool specified_collision_partner = dynamic_cast( + this)->specified_collision_partner; + return std::make_pair(undeclared, specified_collision_partner); + } return std::make_pair(undeclared, false); } @@ -343,9 +386,9 @@ bool Reaction::checkSpecies(const Kinetics& kin) const return true; } -ElementaryReaction::ElementaryReaction(const Composition& reactants_, - const Composition products_, - const Arrhenius& rate_) +ElementaryReaction2::ElementaryReaction2(const Composition& reactants_, + const Composition products_, + const Arrhenius& rate_) : Reaction(reactants_, products_) , rate(rate_) , allow_negative_pre_exponential_factor(false) @@ -353,25 +396,25 @@ ElementaryReaction::ElementaryReaction(const Composition& reactants_, reaction_type = ELEMENTARY_RXN; } -ElementaryReaction::ElementaryReaction() +ElementaryReaction2::ElementaryReaction2() : Reaction() , allow_negative_pre_exponential_factor(false) { reaction_type = ELEMENTARY_RXN; } -void ElementaryReaction::validate() +void ElementaryReaction2::validate() { Reaction::validate(); if (!allow_negative_pre_exponential_factor && rate.preExponentialFactor() < 0) { - throw InputFileError("ElementaryReaction::validate", input, + throw InputFileError("ElementaryReaction2::validate", input, "Undeclared negative pre-exponential factor found in reaction '" + equation() + "'"); } } -void ElementaryReaction::getParameters(AnyMap& reactionNode) const +void ElementaryReaction2::getParameters(AnyMap& reactionNode) const { Reaction::getParameters(reactionNode); if (allow_negative_pre_exponential_factor) { @@ -387,49 +430,62 @@ ThirdBody::ThirdBody(double default_eff) { } +ThirdBody::ThirdBody(const AnyMap& node) +{ + setEfficiencies(node); +} + +void ThirdBody::setEfficiencies(const AnyMap& node) +{ + default_efficiency = node.getDouble("default-efficiency", 1.0); + if (node.hasKey("efficiencies")) { + efficiencies = node["efficiencies"].asMap(); + } +} + double ThirdBody::efficiency(const std::string& k) const { return getValue(efficiencies, k, default_efficiency); } -ThreeBodyReaction::ThreeBodyReaction() +ThreeBodyReaction2::ThreeBodyReaction2() { reaction_type = THREE_BODY_RXN; } -ThreeBodyReaction::ThreeBodyReaction(const Composition& reactants_, - const Composition& products_, - const Arrhenius& rate_, - const ThirdBody& tbody) - : ElementaryReaction(reactants_, products_, rate_) +ThreeBodyReaction2::ThreeBodyReaction2(const Composition& reactants_, + const Composition& products_, + const Arrhenius& rate_, + const ThirdBody& tbody) + : ElementaryReaction2(reactants_, products_, rate_) , third_body(tbody) { reaction_type = THREE_BODY_RXN; } -std::string ThreeBodyReaction::reactantString() const +std::string ThreeBodyReaction2::reactantString() const { if (specified_collision_partner) { - return ElementaryReaction::reactantString() + " + " + return ElementaryReaction2::reactantString() + " + " + third_body.efficiencies.begin()->first; } else { - return ElementaryReaction::reactantString() + " + M"; + return ElementaryReaction2::reactantString() + " + M"; } } -std::string ThreeBodyReaction::productString() const +std::string ThreeBodyReaction2::productString() const { if (specified_collision_partner) { - return ElementaryReaction::productString() + " + " + return ElementaryReaction2::productString() + " + " + third_body.efficiencies.begin()->first; } else { - return ElementaryReaction::productString() + " + M"; + return ElementaryReaction2::productString() + " + M"; } } -void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin) +void ThreeBodyReaction2::calculateRateCoeffUnits(const Kinetics& kin) { - ElementaryReaction::calculateRateCoeffUnits(kin); + ElementaryReaction2::calculateRateCoeffUnits(kin); bool specified_collision_partner_ = false; for (const auto& reac : reactants) { // While this reaction was already identified as a three-body reaction in a @@ -448,9 +504,9 @@ void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin) } } -void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const +void ThreeBodyReaction2::getParameters(AnyMap& reactionNode) const { - ElementaryReaction::getParameters(reactionNode); + ElementaryReaction2::getParameters(reactionNode); if (!specified_collision_partner) { reactionNode["type"] = "three-body"; reactionNode["efficiencies"] = third_body.efficiencies; @@ -461,7 +517,7 @@ void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const } } -std::pair, bool> ThreeBodyReaction::undeclaredThirdBodies( +std::pair, bool> ThreeBodyReaction2::undeclaredThirdBodies( const Kinetics& kin) const { std::vector undeclared; @@ -594,131 +650,47 @@ void ChemicallyActivatedReaction::getParameters(AnyMap& reactionNode) const reactionNode["type"] = "chemically-activated"; } -PlogReaction::PlogReaction() +PlogReaction2::PlogReaction2() : Reaction() { reaction_type = PLOG_RXN; } -PlogReaction::PlogReaction(const Composition& reactants_, - const Composition& products_, const Plog& rate_) +PlogReaction2::PlogReaction2(const Composition& reactants_, + const Composition& products_, const Plog& rate_) : Reaction(reactants_, products_) , rate(rate_) { reaction_type = PLOG_RXN; } -void PlogReaction::getParameters(AnyMap& reactionNode) const +void PlogReaction2::getParameters(AnyMap& reactionNode) const { Reaction::getParameters(reactionNode); reactionNode["type"] = "pressure-dependent-Arrhenius"; - std::vector rateList; - for (const auto& r : rate.rates()) { - AnyMap rateNode; - rateNode["P"].setQuantity(r.first, "Pa"); - r.second.getParameters(rateNode, rate_units); - rateList.push_back(std::move(rateNode)); - } - reactionNode["rate-constants"] = std::move(rateList); + rate.getParameters(reactionNode, rate_units); } -ChebyshevReaction::ChebyshevReaction() +ChebyshevReaction2::ChebyshevReaction2() : Reaction() { reaction_type = CHEBYSHEV_RXN; } -ChebyshevReaction::ChebyshevReaction(const Composition& reactants_, - const Composition& products_, - const ChebyshevRate& rate_) +ChebyshevReaction2::ChebyshevReaction2(const Composition& reactants_, + const Composition& products_, + const Chebyshev& rate_) : Reaction(reactants_, products_) , rate(rate_) { reaction_type = CHEBYSHEV_RXN; } -void Reaction2::setParameters(const AnyMap& node, const Kinetics& kin) -{ - parseReactionEquation(*this, node["equation"], kin); - // Non-stoichiometric reaction orders - std::map orders; - if (node.hasKey("orders")) { - for (const auto& order : node["orders"].asMap()) { - orders[order.first] = order.second; - if (kin.kineticsSpeciesIndex(order.first) == npos) { - setValid(false); - } - } - } - - // Flags - id = node.getString("id", ""); - duplicate = node.getBool("duplicate", false); - allow_negative_orders = node.getBool("negative-orders", false); - allow_nonreactant_orders = node.getBool("nonreactant-orders", false); - - calculateRateCoeffUnits(kin); - input = node; -} - -CustomFunc1Reaction::CustomFunc1Reaction() - : Reaction2() -{ - reaction_type = CUSTOMPY_RXN; -} - -void CustomFunc1Reaction::setParameters(const AnyMap& node, const Kinetics& kin) -{ - Reaction2::setParameters(node, kin); - setRate( - std::shared_ptr(new CustomFunc1Rate(node, rate_units))); -} - -TestReaction::TestReaction() - : Reaction2() - , allow_negative_pre_exponential_factor(false) -{ -} - -void TestReaction::setParameters(const AnyMap& node, const Kinetics& kin) -{ - Reaction2::setParameters(node, kin); - - setRate( - std::shared_ptr(new ArrheniusRate(node, rate_units))); - allow_negative_pre_exponential_factor = node.getBool("negative-A", false); -} - -void ChebyshevReaction::getParameters(AnyMap& reactionNode) const +void ChebyshevReaction2::getParameters(AnyMap& reactionNode) const { Reaction::getParameters(reactionNode); reactionNode["type"] = "Chebyshev"; - reactionNode["temperature-range"].setQuantity({rate.Tmin(), rate.Tmax()}, "K"); - reactionNode["pressure-range"].setQuantity({rate.Pmin(), rate.Pmax()}, "Pa"); - const auto& coeffs1d = rate.coeffs(); - size_t nT = rate.nTemperature(); - size_t nP = rate.nPressure(); - std::vector coeffs2d(nT, vector_fp(nP)); - for (size_t i = 0; i < nT; i++) { - for (size_t j = 0; j < nP; j++) { - coeffs2d[i][j] = coeffs1d[nP*i + j]; - } - } - // Unit conversions must take place later, after the destination unit system - // is known. A lambda function is used here to override the default behavior - Units rate_units2 = rate_units; - auto converter = [rate_units2](AnyValue& coeffs, const UnitSystem& units) { - if (rate_units2.factor() != 0.0) { - coeffs.asVector()[0][0] += std::log10(units.convertFrom(1.0, rate_units2)); - } else if (units.getDelta(UnitSystem()).size()) { - throw CanteraError("ChebyshevReaction::getParameters lambda", - "Cannot convert rate constant with unknown dimensions to a " - "non-default unit system"); - } - }; - AnyValue coeffs; - coeffs = std::move(coeffs2d); - reactionNode["data"].setQuantity(coeffs, converter); + rate.getParameters(reactionNode, rate_units); } InterfaceReaction::InterfaceReaction() @@ -732,7 +704,7 @@ InterfaceReaction::InterfaceReaction(const Composition& reactants_, const Composition& products_, const Arrhenius& rate_, bool isStick) - : ElementaryReaction(reactants_, products_, rate_) + : ElementaryReaction2(reactants_, products_, rate_) , is_sticking_coefficient(isStick) , use_motz_wise_correction(false) { @@ -741,7 +713,7 @@ InterfaceReaction::InterfaceReaction(const Composition& reactants_, void InterfaceReaction::calculateRateCoeffUnits(const Kinetics& kin) { - ElementaryReaction::calculateRateCoeffUnits(kin); + ElementaryReaction2::calculateRateCoeffUnits(kin); if (is_sticking_coefficient || input.hasKey("sticking-coefficient")) { rate_units = Units(1.0); // sticking coefficients are dimensionless } @@ -749,7 +721,7 @@ void InterfaceReaction::calculateRateCoeffUnits(const Kinetics& kin) void InterfaceReaction::getParameters(AnyMap& reactionNode) const { - ElementaryReaction::getParameters(reactionNode); + ElementaryReaction2::getParameters(reactionNode); if (is_sticking_coefficient) { reactionNode["sticking-coefficient"] = std::move(reactionNode["rate-constant"]); reactionNode.erase("rate-constant"); @@ -891,6 +863,242 @@ void BlowersMaselInterfaceReaction::getParameters(AnyMap& reactionNode) const } } +ElementaryReaction3::ElementaryReaction3() +{ + m_rate.reset(new ArrheniusRate); +} + +ElementaryReaction3::ElementaryReaction3(const Composition& reactants, + const Composition& products, + const ArrheniusRate& rate) + : Reaction(reactants, products) +{ + m_rate.reset(new ArrheniusRate(rate)); +} + +ElementaryReaction3::ElementaryReaction3(const AnyMap& node, const Kinetics& kin) + : ElementaryReaction3() +{ + setParameters(node, kin); + setRate(std::make_shared(node, rate_units)); +} + +void ElementaryReaction3::getParameters(AnyMap& reactionNode) const +{ + Reaction::getParameters(reactionNode); + reactionNode.update(m_rate->parameters(rate_units)); +} + +ThreeBodyReaction3::ThreeBodyReaction3() + : ElementaryReaction3() +{ + m_third_body.reset(new ThirdBody); +} + +ThreeBodyReaction3::ThreeBodyReaction3(const Composition& reactants, + const Composition& products, + const ArrheniusRate& rate, + const ThirdBody& tbody) + : ElementaryReaction3(reactants, products, rate) +{ + m_third_body = std::make_shared(tbody); +} + +ThreeBodyReaction3::ThreeBodyReaction3(const AnyMap& node, const Kinetics& kin) + : ThreeBodyReaction3() +{ + setParameters(node, kin); + setRate(std::make_shared(node, rate_units)); +} + +bool ThreeBodyReaction3::detectEfficiencies() +{ + for (const auto& reac : reactants) { + // detect explicitly specified collision partner + if (products.count(reac.first)) { + m_third_body->efficiencies[reac.first] = 1.; + } + } + + if (m_third_body->efficiencies.size() == 0) { + return false; + } else if (m_third_body->efficiencies.size() > 1) { + throw CanteraError("ThreeBodyReaction3::detectEfficiencies", + "Found more than one explicitly specified collision partner\n" + "in reaction '{}'.", equation()); + } + + m_third_body->default_efficiency = 0.; + 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 ThreeBodyReaction3::calculateRateCoeffUnits(const Kinetics& kin) +{ + ElementaryReaction3::calculateRateCoeffUnits(kin); + bool specified_collision_partner_ = false; + for (const auto& reac : reactants) { + // While this reaction was already identified as a three-body reaction in a + // pre-processing step, this method is often called before a three-body + // reaction is fully instantiated. For the determination of the correct units, + // it is necessary to check whether the reaction uses a generic 'M' or an + // explicitly specified collision partner that may not have been deleted yet. + if (reac.first != "M" && products.count(reac.first)) { + // detected specified third-body collision partner + specified_collision_partner_ = true; + } + } + if (!specified_collision_partner_) { + const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex()); + rate_units *= rxn_phase.standardConcentrationUnits().pow(-1); + } +} + +void ThreeBodyReaction3::setParameters(const AnyMap& node, const Kinetics& kin) +{ + if (node.empty()) { + // empty node: used by newReaction() factory loader + return; + } + Reaction::setParameters(node, kin); + if (reactants.count("M") != 1 || products.count("M") != 1) { + if (!detectEfficiencies()) { + throw InputFileError("ThreeBodyReaction3::setParameters", node["equation"], + "Reaction equation '{}' does not contain third body 'M'", + node["equation"].asString()); + } + return; + } + + reactants.erase("M"); + products.erase("M"); + m_third_body->setEfficiencies(node); +} + +void ThreeBodyReaction3::getParameters(AnyMap& reactionNode) const +{ + ElementaryReaction3::getParameters(reactionNode); + if (!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::string ThreeBodyReaction3::reactantString() const +{ + if (specified_collision_partner) { + return ElementaryReaction3::reactantString() + " + " + + m_third_body->efficiencies.begin()->first; + } else { + return ElementaryReaction3::reactantString() + " + M"; + } +} + +std::string ThreeBodyReaction3::productString() const +{ + if (specified_collision_partner) { + return ElementaryReaction3::productString() + " + " + + m_third_body->efficiencies.begin()->first; + } else { + return ElementaryReaction3::productString() + " + M"; + } +} + +PlogReaction3::PlogReaction3() +{ + m_rate.reset(new PlogRate); +} + +PlogReaction3::PlogReaction3(const Composition& reactants, + const Composition& products, const PlogRate& rate) + : Reaction(reactants, products) +{ + m_rate.reset(new PlogRate(rate)); +} + +PlogReaction3::PlogReaction3(const AnyMap& node, const Kinetics& kin) + : PlogReaction3() +{ + setParameters(node, kin); + setRate(std::make_shared(node, rate_units)); +} + +void PlogReaction3::getParameters(AnyMap& reactionNode) const +{ + Reaction::getParameters(reactionNode); + reactionNode["type"] = "pressure-dependent-Arrhenius"; + reactionNode.update(m_rate->parameters(rate_units)); +} + +ChebyshevReaction3::ChebyshevReaction3() +{ + m_rate.reset(new ChebyshevRate3); +} + +ChebyshevReaction3::ChebyshevReaction3(const Composition& reactants, + const Composition& products, + const ChebyshevRate3& rate) + : Reaction(reactants, products) +{ + m_rate.reset(new ChebyshevRate3(rate)); +} + +ChebyshevReaction3::ChebyshevReaction3(const AnyMap& node, const Kinetics& kin) + : ChebyshevReaction3() +{ + setParameters(node, kin); + setRate(std::make_shared(node, rate_units)); +} + +void ChebyshevReaction3::getParameters(AnyMap& reactionNode) const +{ + Reaction::getParameters(reactionNode); + reactionNode["type"] = "Chebyshev"; + reactionNode.update(m_rate->parameters(rate_units)); +} + +CustomFunc1Reaction::CustomFunc1Reaction() +{ + m_rate.reset(new CustomFunc1Rate); +} + +CustomFunc1Reaction::CustomFunc1Reaction(const Composition& reactants, + const Composition& products, + const CustomFunc1Rate& rate) + : Reaction(reactants, products) +{ + m_rate.reset(new CustomFunc1Rate(rate)); +} + +CustomFunc1Reaction::CustomFunc1Reaction(const AnyMap& node, const Kinetics& kin) + : CustomFunc1Reaction() +{ + setParameters(node, kin); + setRate(std::make_shared(node, rate_units)); +} + Arrhenius readArrhenius(const XML_Node& arrhenius_node) { return Arrhenius(getFloat(arrhenius_node, "A", "toSI"), @@ -1008,10 +1216,7 @@ void readEfficiencies(ThirdBody& tbody, const XML_Node& rc_node) void readEfficiencies(ThirdBody& tbody, const AnyMap& node) { - tbody.default_efficiency = node.getDouble("default-efficiency", 1.0); - if (node.hasKey("efficiencies")) { - tbody.efficiencies = node["efficiencies"].asMap(); - } + tbody.setEfficiencies(node); } BlowersMasel readBlowersMasel(const Reaction& R, const AnyValue& rate, @@ -1040,7 +1245,7 @@ BlowersMasel readBlowersMasel(const Reaction& R, const AnyValue& rate, return BlowersMasel(A, b, Ta0, w); } -bool detectEfficiencies(ThreeBodyReaction& R) +bool detectEfficiencies(ThreeBodyReaction2& R) { for (const auto& reac : R.reactants) { // detect explicitly specified collision partner @@ -1168,7 +1373,6 @@ void setupReaction(Reaction& R, const AnyMap& node, const Kinetics& kin) { parseReactionEquation(R, node["equation"], kin); // Non-stoichiometric reaction orders - std::map orders; if (node.hasKey("orders")) { for (const auto& order : node["orders"].asMap()) { R.orders[order.first] = order.second; @@ -1188,7 +1392,7 @@ void setupReaction(Reaction& R, const AnyMap& node, const Kinetics& kin) R.calculateRateCoeffUnits(kin); } -void setupElementaryReaction(ElementaryReaction& R, const XML_Node& rxn_node) +void setupElementaryReaction(ElementaryReaction2& R, const XML_Node& rxn_node) { const XML_Node& rc_node = rxn_node.child("rateCoeff"); if (rc_node.hasChild("Arrhenius")) { @@ -1210,7 +1414,7 @@ void setupElementaryReaction(ElementaryReaction& R, const XML_Node& rxn_node) setupReaction(R, rxn_node); } -void setupElementaryReaction(ElementaryReaction& R, const AnyMap& node, +void setupElementaryReaction(ElementaryReaction2& R, const AnyMap& node, const Kinetics& kin) { setupReaction(R, node, kin); @@ -1218,7 +1422,7 @@ void setupElementaryReaction(ElementaryReaction& R, const AnyMap& node, R.rate = readArrhenius(R, node["rate-constant"], kin, node.units()); } -void setupThreeBodyReaction(ThreeBodyReaction& R, const XML_Node& rxn_node) +void setupThreeBodyReaction(ThreeBodyReaction2& R, const XML_Node& rxn_node) { readEfficiencies(R.third_body, rxn_node.child("rateCoeff")); setupElementaryReaction(R, rxn_node); @@ -1227,7 +1431,7 @@ void setupThreeBodyReaction(ThreeBodyReaction& R, const XML_Node& rxn_node) } } -void setupThreeBodyReaction(ThreeBodyReaction& R, const AnyMap& node, +void setupThreeBodyReaction(ThreeBodyReaction2& R, const AnyMap& node, const Kinetics& kin) { setupElementaryReaction(R, node, kin); @@ -1350,7 +1554,7 @@ void setupChemicallyActivatedReaction(ChemicallyActivatedReaction& R, setupReaction(R, rxn_node); } -void setupPlogReaction(PlogReaction& R, const XML_Node& rxn_node) +void setupPlogReaction(PlogReaction2& R, const XML_Node& rxn_node) { XML_Node& rc = rxn_node.child("rateCoeff"); std::multimap rates; @@ -1362,7 +1566,7 @@ void setupPlogReaction(PlogReaction& R, const XML_Node& rxn_node) setupReaction(R, rxn_node); } -void setupPlogReaction(PlogReaction& R, const AnyMap& node, const Kinetics& kin) +void setupPlogReaction(PlogReaction2& R, const AnyMap& node, const Kinetics& kin) { setupReaction(R, node, kin); std::multimap rates; @@ -1373,13 +1577,13 @@ void setupPlogReaction(PlogReaction& R, const AnyMap& node, const Kinetics& kin) R.rate = Plog(rates); } -void PlogReaction::validate() +void PlogReaction2::validate() { Reaction::validate(); rate.validate(equation()); } -void setupChebyshevReaction(ChebyshevReaction& R, const XML_Node& rxn_node) +void setupChebyshevReaction(ChebyshevReaction2& R, const XML_Node& rxn_node) { XML_Node& rc = rxn_node.child("rateCoeff"); const XML_Node& coeff_node = rc.child("floatArray"); @@ -1394,15 +1598,15 @@ void setupChebyshevReaction(ChebyshevReaction& R, const XML_Node& rxn_node) coeffs(t,p) = coeffs_flat[nP*t + p]; } } - R.rate = ChebyshevRate(getFloat(rc, "Tmin", "toSI"), - getFloat(rc, "Tmax", "toSI"), - getFloat(rc, "Pmin", "toSI"), - getFloat(rc, "Pmax", "toSI"), - coeffs); + R.rate = Chebyshev(getFloat(rc, "Tmin", "toSI"), + getFloat(rc, "Tmax", "toSI"), + getFloat(rc, "Pmin", "toSI"), + getFloat(rc, "Pmax", "toSI"), + coeffs); setupReaction(R, rxn_node); } -void setupChebyshevReaction(ChebyshevReaction&R, const AnyMap& node, +void setupChebyshevReaction(ChebyshevReaction2&R, const AnyMap& node, const Kinetics& kin) { setupReaction(R, node, kin); @@ -1423,11 +1627,11 @@ void setupChebyshevReaction(ChebyshevReaction&R, const AnyMap& node, } const UnitSystem& units = node.units(); coeffs(0, 0) += std::log10(units.convertTo(1.0, R.rate_units)); - R.rate = ChebyshevRate(units.convert(T_range[0], "K"), - units.convert(T_range[1], "K"), - units.convert(P_range[0], "Pa"), - units.convert(P_range[1], "Pa"), - coeffs); + R.rate = Chebyshev(units.convert(T_range[0], "K"), + units.convert(T_range[1], "K"), + units.convert(P_range[0], "Pa"), + units.convert(P_range[1], "Pa"), + coeffs); } void setupInterfaceReaction(InterfaceReaction& R, const XML_Node& rxn_node) diff --git a/src/kinetics/ReactionData.cpp b/src/kinetics/ReactionData.cpp index 07041e33a21..9da9832f91c 100644 --- a/src/kinetics/ReactionData.cpp +++ b/src/kinetics/ReactionData.cpp @@ -5,15 +5,40 @@ #include "cantera/kinetics/ReactionData.h" #include "cantera/thermo/ThermoPhase.h" +#include "cantera/base/ctexceptions.h" namespace Cantera { -void ArrheniusData::update(const ThermoPhase& bulk) { +void ArrheniusData::update(const ThermoPhase& bulk) +{ update(bulk.temperature()); } -void CustomFunc1Data::update(const ThermoPhase& bulk) { +void PlogData::update(double T) +{ + throw CanteraError("PlogData::update", + "Missing state information: reaction type requires pressure."); +} + +void PlogData::update(const ThermoPhase& bulk) +{ + update(bulk.temperature(), bulk.pressure()); +} + +void ChebyshevData::update(double T) +{ + throw CanteraError("ChebyshevData::update", + "Missing state information: reaction type requires pressure."); +} + +void ChebyshevData::update(const ThermoPhase& bulk) +{ + update(bulk.temperature(), bulk.pressure()); +} + +void CustomFunc1Data::update(const ThermoPhase& bulk) +{ m_temperature = bulk.temperature(); } diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index b0b7834ac69..bdffc49332a 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -22,141 +22,204 @@ std::mutex ReactionFactory::reaction_mutex; ReactionFactory::ReactionFactory() { // register elementary reactions - reg("elementary", []() { return new ElementaryReaction(); }); + reg("elementary", [](const AnyMap& node, const Kinetics& kin) { + return new ElementaryReaction3(node, kin); + }); addAlias("elementary", "arrhenius"); addAlias("elementary", ""); - reg_XML("elementary", - [](Reaction* R, const XML_Node& node) { - setupElementaryReaction(*(ElementaryReaction*)R, node); - }); - reg_AnyMap("elementary", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupElementaryReaction(*(ElementaryReaction*)R, node, kin); - }); + + // register elementary reactions (old framework) + reg("elementary-legacy", [](const AnyMap& node, const Kinetics& kin) { + ElementaryReaction* R = new ElementaryReaction2(); + if (!node.empty()) { + setupElementaryReaction(*R, node, kin); + } + return R; + }); // register three-body reactions - reg("three-body", []() { return new ThreeBodyReaction(); }); + reg("three-body", [](const AnyMap& node, const Kinetics& kin) { + return new ThreeBodyReaction3(node, kin); + }); addAlias("three-body", "threebody"); addAlias("three-body", "three_body"); - reg_XML("three-body", - [](Reaction* R, const XML_Node& node) { - setupThreeBodyReaction(*(ThreeBodyReaction*)R, node); - }); - reg_AnyMap("three-body", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupThreeBodyReaction(*(ThreeBodyReaction*)R, node, kin); - }); + + // register three-body reactions (old framework) + reg("three-body-legacy", [](const AnyMap& node, const Kinetics& kin) { + ThreeBodyReaction* R = new ThreeBodyReaction2(); + if (!node.empty()) { + setupThreeBodyReaction(*R, node, kin); + } + return R; + }); // register falloff reactions - reg("falloff", []() { return new FalloffReaction(); }); - reg_XML("falloff", - [](Reaction* R, const XML_Node& node) { - setupFalloffReaction(*(FalloffReaction*)R, node); - }); - reg_AnyMap("falloff", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupFalloffReaction(*(FalloffReaction*)R, node, kin); - }); + reg("falloff", [](const AnyMap& node, const Kinetics& kin) { + FalloffReaction* R = new FalloffReaction(); + if (!node.empty()) { + setupFalloffReaction(*R, node, kin); + } + return R; + }); // register falloff reactions - reg("chemically-activated", []() { return new ChemicallyActivatedReaction(); }); + reg("chemically-activated", [](const AnyMap& node, const Kinetics& kin) { + FalloffReaction* R = new ChemicallyActivatedReaction(); + if (!node.empty()) { + setupFalloffReaction(*R, node, kin); + } + return R; + }); addAlias("chemically-activated", "chemact"); addAlias("chemically-activated", "chemically_activated"); - reg_XML("chemically-activated", - [](Reaction* R, const XML_Node& node) { - setupChemicallyActivatedReaction(*(ChemicallyActivatedReaction*)R, node); - }); - reg_AnyMap("chemically-activated", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupFalloffReaction(*(FalloffReaction*)R, node, kin); - }); - // register pressure-depdendent-Arrhenius reactions - reg("pressure-dependent-Arrhenius", []() { return new PlogReaction(); }); + // register pressure-dependent-Arrhenius reactions + reg("pressure-dependent-Arrhenius", [](const AnyMap& node, const Kinetics& kin) { + return new PlogReaction3(node, kin); + }); addAlias("pressure-dependent-Arrhenius", "plog"); addAlias("pressure-dependent-Arrhenius", "pdep_arrhenius"); - reg_XML("pressure-dependent-Arrhenius", - [](Reaction* R, const XML_Node& node) { - setupPlogReaction(*(PlogReaction*)R, node); - }); - reg_AnyMap("pressure-dependent-Arrhenius", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupPlogReaction(*(PlogReaction*)R, node, kin); - }); + + // register pressure-dependent-Arrhenius reactions (old framework) + reg("pressure-dependent-Arrhenius-legacy", [](const AnyMap& node, const Kinetics& kin) { + PlogReaction* R = new PlogReaction2(); + if (!node.empty()) { + setupPlogReaction(*R, node, kin); + } + return R; + }); // register Chebyshev reactions - reg("Chebyshev", []() { return new ChebyshevReaction(); }); + reg("Chebyshev", [](const AnyMap& node, const Kinetics& kin) { + return new ChebyshevReaction3(node, kin); + }); addAlias("Chebyshev", "chebyshev"); - reg_XML("Chebyshev", - [](Reaction* R, const XML_Node& node) { - setupChebyshevReaction(*(ChebyshevReaction*)R, node); - }); - reg_AnyMap("Chebyshev", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupChebyshevReaction(*(ChebyshevReaction*)R, node, kin); - }); + reg("Chebyshev-legacy", [](const AnyMap& node, const Kinetics& kin) { + ChebyshevReaction* R = new ChebyshevReaction2(); + if (!node.empty()) { + setupChebyshevReaction(*R, node, kin); + } + return R; + }); // register custom reactions specified by Func1 objects - reg("custom-rate-function", []() { return new CustomFunc1Reaction(); }); - reg_XML("custom-rate-function", - [](Reaction* R, const XML_Node& node) { - throw CanteraError("ReactionFactory", "Custom reactions based " - "on 'Func1' objects cannot be created from XML nodes'"); - }); - reg_AnyMap("custom-rate-function", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - ((Reaction2*)R)->setParameters(node, kin); - }); - - // register custom Python reactions - reg("elementary-new", []() { return new TestReaction(); }); - reg_XML("elementary-new", - [](Reaction* R, const XML_Node& node) { - throw CanteraError("ReactionFactory", "Test reactions " - "cannot be created from XML nodes'"); - }); - reg_AnyMap("elementary-new", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - (dynamic_cast(R))->setParameters(node, kin); - }); + reg("custom-rate-function", [](const AnyMap& node, const Kinetics& kin) { + return new CustomFunc1Reaction(node, kin); + }); // register interface reactions - reg("interface", []() { return new InterfaceReaction(); }); + reg("interface", [](const AnyMap& node, const Kinetics& kin) { + InterfaceReaction* R = new InterfaceReaction(); + if (!node.empty()) { + setupInterfaceReaction(*R, node, kin); + } + return R; + }); addAlias("interface", "surface"); addAlias("interface", "edge"); - reg_XML("interface", - [](Reaction* R, const XML_Node& node) { - setupInterfaceReaction(*(InterfaceReaction*)R, node); - }); - reg_AnyMap("interface", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupInterfaceReaction(*(InterfaceReaction*)R, node, kin); - }); // register electrochemical reactions - reg("electrochemical", []() { return new ElectrochemicalReaction(); }); - reg_XML("electrochemical", - [](Reaction* R, const XML_Node& node) { - setupElectrochemicalReaction(*(ElectrochemicalReaction*)R, node); - }); - reg_AnyMap("electrochemical", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupElectrochemicalReaction(*(ElectrochemicalReaction*)R, node, kin); - }); + reg("electrochemical", [](const AnyMap& node, const Kinetics& kin) { + ElectrochemicalReaction* R = new ElectrochemicalReaction(); + if (!node.empty()) { + setupElectrochemicalReaction(*R, node, kin); + } + return R; + }); // register Blowers Masel reactions - reg("Blowers-Masel", []() { return new BlowersMaselReaction(); }); - reg_AnyMap("Blowers-Masel", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupBlowersMaselReaction(*(BlowersMaselReaction*)R, node, kin); - }); + reg("Blowers-Masel", [](const AnyMap& node, const Kinetics& kin) { + BlowersMaselReaction* R = new BlowersMaselReaction(); + if (!node.empty()) { + setupBlowersMaselReaction(*R, node, kin); + } + return R; + }); // register surface Blowers Masel reactions - reg("surface-Blowers-Masel", []() { return new BlowersMaselInterfaceReaction(); }); - reg_AnyMap("surface-Blowers-Masel", - [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupBlowersMaselInterfaceReaction(*(BlowersMaselInterfaceReaction*)R, node, kin); - }); + reg("surface-Blowers-Masel", [](const AnyMap& node, const Kinetics& kin) { + BlowersMaselInterfaceReaction* R = new BlowersMaselInterfaceReaction(); + if (!node.empty()) { + setupBlowersMaselInterfaceReaction(*R, node, kin); + } + return R; + }); +} + +ReactionFactoryXML* ReactionFactoryXML::s_factory = 0; +std::mutex ReactionFactoryXML::reaction_mutex; + +ReactionFactoryXML::ReactionFactoryXML() +{ + // register elementary reactions + reg("elementary-legacy", [](const XML_Node& node) { + Reaction* R = new ElementaryReaction2(); + setupElementaryReaction(*(ElementaryReaction2*)R, node); + return R; + }); + addAlias("elementary-legacy", "elementary"); + addAlias("elementary-legacy", "arrhenius"); + addAlias("elementary-legacy", ""); + + // register three-body reactions + reg("three-body-legacy", [](const XML_Node& node) { + Reaction* R = new ThreeBodyReaction2(); + setupThreeBodyReaction(*(ThreeBodyReaction2*)R, node); + return R; + }); + addAlias("three-body-legacy", "three-body"); + addAlias("three-body-legacy", "threebody"); + addAlias("three-body-legacy", "three_body"); + + // register falloff reactions + reg("falloff", [](const XML_Node& node) { + Reaction* R = new FalloffReaction(); + setupFalloffReaction(*(FalloffReaction*)R, node); + return R; + }); + + // register falloff reactions + reg("chemically-activated", [](const XML_Node& node) { + Reaction* R = new ChemicallyActivatedReaction(); + setupChemicallyActivatedReaction(*(ChemicallyActivatedReaction*)R, node); + return R; + }); + addAlias("chemically-activated", "chemact"); + addAlias("chemically-activated", "chemically_activated"); + + // register pressure-depdendent-Arrhenius reactions + reg("pressure-dependent-Arrhenius-legacy", [](const XML_Node& node) { + Reaction* R = new PlogReaction2(); + setupPlogReaction(*(PlogReaction2*)R, node); + return R; + }); + addAlias("pressure-dependent-Arrhenius-legacy", "pressure-dependent-Arrhenius"); + addAlias("pressure-dependent-Arrhenius-legacy", "plog"); + addAlias("pressure-dependent-Arrhenius-legacy", "pdep_arrhenius"); + + // register Chebyshev reactions + reg("Chebyshev-legacy", [](const XML_Node& node) { + Reaction* R = new ChebyshevReaction2(); + setupChebyshevReaction(*(ChebyshevReaction2*)R, node); + return R; + }); + addAlias("Chebyshev-legacy", "chebyshev"); + + // register interface reactions + reg("interface", [](const XML_Node& node) { + Reaction* R = new InterfaceReaction(); + setupInterfaceReaction(*(InterfaceReaction*)R, node); + return R; + }); + addAlias("interface", "surface"); + addAlias("interface", "edge"); + + // register electrochemical reactions + reg("electrochemical", [](const XML_Node& node) { + Reaction* R = new ElectrochemicalReaction(); + setupElectrochemicalReaction(*(ElectrochemicalReaction*)R, node); + return R; + }); } bool isThreeBody(const Reaction& R) @@ -228,7 +291,9 @@ bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin) unique_ptr newReaction(const std::string& type) { - unique_ptr R(ReactionFactory::factory()->create(type)); + AnyMap rxn_node; + Kinetics kin; + unique_ptr R(ReactionFactory::factory()->create(type, rxn_node, kin)); return R; } @@ -243,33 +308,24 @@ unique_ptr newReaction(const XML_Node& rxn_node) type = "electrochemical"; } - Reaction* R; - try { - R = ReactionFactory::factory()->create(type); - } catch (CanteraError& err) { + if (!(ReactionFactoryXML::factory()->exists(type))) { throw CanteraError("newReaction", "Unknown reaction type '" + rxn_node["type"] + "'"); } if (type.empty()) { // Reaction type is not specified // See if this is a three-body reaction with a specified collision partner - ElementaryReaction testReaction; + ElementaryReaction2 testReaction; setupReaction(testReaction, rxn_node); if (isThreeBody(testReaction)) { type = "three-body"; - R = ReactionFactory::factory()->create(type); } } - if (type != "electrochemical") { - type = R->type(); - } - ReactionFactory::factory()->setup_XML(type, R, rxn_node); - + Reaction* R = ReactionFactoryXML::factory()->create(type, rxn_node); return unique_ptr(R); } -unique_ptr newReaction(const AnyMap& rxn_node, - const Kinetics& kin) +unique_ptr newReaction(const AnyMap& rxn_node, const Kinetics& kin) { std::string type = "elementary"; if (rxn_node.hasKey("type")) { @@ -277,7 +333,7 @@ unique_ptr newReaction(const AnyMap& rxn_node, } else if (kin.thermo().nDim() == 3) { // Reaction type is not specified // See if this is a three-body reaction with a specified collision partner - ElementaryReaction testReaction; + ElementaryReaction2 testReaction; parseReactionEquation(testReaction, rxn_node["equation"], kin); if (isThreeBody(testReaction)) { type = "three-body"; @@ -287,7 +343,7 @@ unique_ptr newReaction(const AnyMap& rxn_node, if (kin.thermo().nDim() < 3 && type == "elementary") { // See if this is an electrochemical reaction: type of // receiving reaction object is unimportant in this case - ElementaryReaction testReaction; + ElementaryReaction2 testReaction; parseReactionEquation(testReaction, rxn_node["equation"], kin); if (isElectrochemicalReaction(testReaction, kin)) { type = "electrochemical"; @@ -300,15 +356,11 @@ unique_ptr newReaction(const AnyMap& rxn_node, type = "surface-Blowers-Masel"; } - Reaction* R; - try { - R = ReactionFactory::factory()->create(type); - } catch (CanteraError& err) { + if (!(ReactionFactory::factory()->exists(type))) { throw InputFileError("ReactionFactory::newReaction", rxn_node["type"], "Unknown reaction type '{}'", type); } - ReactionFactory::factory()->setup_AnyMap(type, R, rxn_node, kin); - + Reaction* R = ReactionFactory::factory()->create(type, rxn_node, kin); return unique_ptr(R); } diff --git a/src/kinetics/ReactionRate.cpp b/src/kinetics/ReactionRate.cpp index 29485704c97..695d210a84e 100644 --- a/src/kinetics/ReactionRate.cpp +++ b/src/kinetics/ReactionRate.cpp @@ -5,31 +5,183 @@ #include "cantera/kinetics/ReactionRate.h" #include "cantera/numerics/Func1.h" -#include "cantera/base/AnyMap.h" namespace Cantera { +void ReactionRateBase::setParameters(const AnyMap& node, const Units& rate_units) +{ + units = rate_units; + input = node; +} + +AnyMap ReactionRateBase::parameters(const Units& rate_units) const +{ + AnyMap out; + getParameters(out, rate_units); + return out; +} + +AnyMap ReactionRateBase::parameters() const +{ + AnyMap out; + getParameters(out, units); + return out; +} + +ArrheniusRate::ArrheniusRate() + : Arrhenius(NAN, NAN, NAN) + , allow_negative_pre_exponential_factor(false) +{ +} + ArrheniusRate::ArrheniusRate(double A, double b, double E) - : Arrhenius(A, b, E) { + : Arrhenius(A, b, E / GasConstant) + , allow_negative_pre_exponential_factor(false) +{ +} + +ArrheniusRate::ArrheniusRate(const AnyMap& node, const Units& rate_units) +{ + setParameters(node, rate_units); +} + +ArrheniusRate::ArrheniusRate(const AnyMap& node) +{ + setParameters(node, Units(1.)); +} + +ArrheniusRate::ArrheniusRate(const Arrhenius& arr, bool allow_negative_A) + : Arrhenius(arr.preExponentialFactor(), + arr.temperatureExponent(), + arr.activationEnergy_R()) + , allow_negative_pre_exponential_factor(allow_negative_A) +{ +} + +void ArrheniusRate::setParameters(const AnyMap& node, const Units& rate_units) +{ + ReactionRateBase::setParameters(node, rate_units); + allow_negative_pre_exponential_factor = node.getBool("negative-A", false); + if (!node.hasKey("rate-constant")) { + Arrhenius::setParameters(AnyValue(), node.units(), rate_units); + return; + } + Arrhenius::setParameters(node["rate-constant"], node.units(), rate_units); +} + +void ArrheniusRate::getParameters(AnyMap& rateNode, + const Units& rate_units) const +{ + if (allow_negative_pre_exponential_factor) { + rateNode["negative-A"] = true; + } + AnyMap node; + Arrhenius::getParameters(node, rate_units); + if (!node.empty()) { + // Arrhenius object is configured + rateNode["rate-constant"] = std::move(node); + } +} + +void ArrheniusRate::validate(const std::string& equation) +{ + if (!allow_negative_pre_exponential_factor && preExponentialFactor() < 0) { + throw CanteraError("ArrheniusRate::validate", + "Undeclared negative pre-exponential factor found in reaction '" + + equation + "'"); + } } -ArrheniusRate::ArrheniusRate(const AnyMap& node, const Units& rate_units) { - setParameters(node["rate-constant"], node.units(), rate_units); +PlogRate::PlogRate(const std::multimap& rates) + : Plog(rates) +{ +} + +PlogRate::PlogRate(const AnyMap& node, const Units& rate_units) +{ + setParameters(node, rate_units); +} + +PlogRate::PlogRate(const AnyMap& node) +{ + setParameters(node, Units(1.)); +} + +void PlogRate::setParameters(const AnyMap& node, const Units& rate_units) +{ + // @TODO implementation of Plog::setParameters should be transferred here + // when the Plog class is removed from RxnRates.h after Cantera 2.6 + ReactionRateBase::setParameters(node, rate_units); + if (!node.hasKey("rate-constants")) { + Plog::setParameters(std::vector (), node.units(), rate_units); + return; + } + Plog::setParameters(node.at("rate-constants").asVector(), + node.units(), rate_units); +} + +void PlogRate::getParameters(AnyMap& rateNode, const Units& rate_units) const +{ + // @TODO implementation of Plog::getParameters should be transferred here + // when the Plog class is removed from RxnRates.h after Cantera 2.6 + Plog::getParameters(rateNode, rate_units); +} + +ChebyshevRate3::ChebyshevRate3(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs) + : Chebyshev(Tmin, Tmax, Pmin, Pmax, coeffs) +{ +} + +ChebyshevRate3::ChebyshevRate3(const AnyMap& node, const Units& rate_units) +{ + setParameters(node, rate_units); +} + +ChebyshevRate3::ChebyshevRate3(const AnyMap& node) +{ + setParameters(node, Units(1.)); +} + +void ChebyshevRate3::setParameters(const AnyMap& node, const Units& rate_units) +{ + ReactionRateBase::setParameters(node, rate_units); + if (!node.hasKey("data")) { + Chebyshev::setParameters(AnyMap(), node.units(), rate_units); + return; + } + // @TODO implementation of Chebyshev::setParameters should be transferred here + // when the Chebyshev class is removed from RxnRates.h after Cantera 2.6 + Chebyshev::setParameters(node, node.units(), rate_units); +} + +void ChebyshevRate3::getParameters(AnyMap& rateNode, + const Units& rate_units) const +{ + // @TODO implementation of Chebyshev::getParameters should be transferred here + // when the Chebyshev class is removed from RxnRates.h after Cantera 2.6 + Chebyshev::getParameters(rateNode, rate_units); +} + +void ChebyshevRate3::validate(const std::string& equation) +{ } CustomFunc1Rate::CustomFunc1Rate() : m_ratefunc(0) {} -void CustomFunc1Rate::setRateFunction(shared_ptr f) { +void CustomFunc1Rate::setRateFunction(shared_ptr f) +{ m_ratefunc = f; } -double CustomFunc1Rate::eval(const CustomFunc1Data& shared_data) const { +double CustomFunc1Rate::eval(const CustomFunc1Data& shared_data, + double concm) const +{ if (m_ratefunc) { return m_ratefunc->eval(shared_data.m_temperature); } - throw CanteraError("CustomFunc1Rate::eval", - "Custom rate function is not initialized."); + return NAN; } } diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 5571aca3042..892821dc8ec 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -6,6 +6,7 @@ #include "cantera/kinetics/RxnRates.h" #include "cantera/base/Array.h" #include "cantera/base/AnyMap.h" +#include "cantera/base/global.h" namespace Cantera { @@ -29,10 +30,20 @@ Arrhenius::Arrhenius(doublereal A, doublereal b, doublereal E) } } +Arrhenius::Arrhenius(const AnyValue& rate, + const UnitSystem& units, const Units& rate_units) +{ + setParameters(rate, units, rate_units); +} + void Arrhenius::setParameters(const AnyValue& rate, const UnitSystem& units, const Units& rate_units) { - if (rate.is()) { + if (rate.empty()) { + m_A = NAN; + m_b = NAN; + m_E = NAN; + } else if (rate.is()) { auto& rate_map = rate.as(); m_A = units.convert(rate_map["A"], rate_units); m_b = rate_map["b"].asDouble(); @@ -44,7 +55,7 @@ void Arrhenius::setParameters(const AnyValue& rate, m_E = units.convertActivationEnergy(rate_vec[2], "K"); } - if (m_A <= 0.0) { + if (m_A <= 0.0) { m_logA = -1.0E300; } else { m_logA = std::log(m_A); @@ -53,10 +64,14 @@ void Arrhenius::setParameters(const AnyValue& rate, void Arrhenius::getParameters(AnyMap& rateNode, const Units& rate_units) const { - if (rate_units.factor() != 0.0) { - rateNode["A"].setQuantity(preExponentialFactor(), rate_units); + double A = preExponentialFactor(); + if (std::isnan(A)) { + // Return empty/unmodified AnyMap + return; + } else if (rate_units.factor() != 0.0) { + rateNode["A"].setQuantity(A, rate_units); } else { - rateNode["A"] = preExponentialFactor(); + rateNode["A"] = A; // This can't be converted to a different unit system because the dimensions of // the rate constant were not set. Can occur if the reaction was created outside // the context of a Kinetics object and never added to a Kinetics object. @@ -140,11 +155,55 @@ void SurfaceArrhenius::addCoverageDependence(size_t k, doublereal a, } } -Plog::Plog(const std::multimap& rates) +Plog::Plog() : logP_(-1000) , logP1_(1000) , logP2_(-1000) , rDeltaP_(-1.0) +{ +} + +Plog::Plog(const std::multimap& rates) + : Plog() +{ + setup(rates); +} + +void Plog::setParameters(const std::vector& rates, + const UnitSystem& units, const Units& rate_units) +{ + std::multimap multi_rates; + if (rates.size()) { + for (const auto& rate : rates) { + multi_rates.insert({rate.convert("P", "Pa"), + Arrhenius(AnyValue(rate), units, rate_units)}); + } + } else { + // ensure that reaction rate can be evaluated (but returns NaN) + multi_rates.insert({1.e-7, Arrhenius(NAN, NAN, NAN)}); + multi_rates.insert({1.e14, Arrhenius(NAN, NAN, NAN)}); + } + setup(multi_rates); +} + +void Plog::getParameters(AnyMap& rateNode, const Units& rate_units) const +{ + std::vector rateList; + double A = rates_[1].preExponentialFactor(); + if (std::isnan(A)) { + // Return empty/unmodified AnyMap + return; + } + for (const auto& r : rates()) { + AnyMap rateNode_; + rateNode_["P"].setQuantity(r.first, "Pa"); + r.second.getParameters(rateNode_, rate_units); + rateList.push_back(std::move(rateNode_)); + } + rateNode["rate-constants"] = std::move(rateList); +} + +void Plog::setup(const std::multimap& rates) { size_t j = 0; rates_.reserve(rates.size()); @@ -178,7 +237,7 @@ void Plog::validate(const std::string& equation) for (size_t i=0; i < 6; i++) { double k = 0; for (size_t p = ilow1_; p < ilow2_; p++) { - k += rates_[p].updateRC(log(T[i]), 1.0/T[i]); + k += rates_.at(p).updateRC(log(T[i]), 1.0/T[i]); } if (k < 0) { format_to(err_reactions, @@ -209,9 +268,8 @@ std::vector > Plog::rates() const return R; } - -ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, - const Array2D& coeffs) +Chebyshev::Chebyshev(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs) : Tmin_(Tmin) , Tmax_(Tmax) , Pmin_(Pmin) @@ -220,6 +278,55 @@ ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, , nT_(coeffs.nRows()) , chebCoeffs_(coeffs.nColumns() * coeffs.nRows(), 0.0) , dotProd_(coeffs.nRows()) +{ + setup(Tmin, Tmax, Pmin, Pmax, coeffs); +} + +void Chebyshev::setParameters(const AnyMap& node, + const UnitSystem& units, const Units& rate_units) +{ + Array2D coeffs; + if (!node.empty()) { + const auto& T_range = node["temperature-range"].asVector(2); + const auto& P_range = node["pressure-range"].asVector(2); + auto& vcoeffs = node["data"].asVector(); + coeffs = Array2D(vcoeffs.size(), vcoeffs[0].size()); + for (size_t i = 0; i < coeffs.nRows(); i++) { + if (vcoeffs[i].size() != vcoeffs[0].size()) { + throw InputFileError("Chebyshev::setParameters", node["data"], + "Inconsistent number of coefficients in row {} of matrix", i + 1); + } + for (size_t j = 0; j < coeffs.nColumns(); j++) { + coeffs(i, j) = vcoeffs[i][j]; + } + } + coeffs(0, 0) += std::log10(units.convertTo(1.0, rate_units)); + + Tmin_ = units.convert(T_range[0], "K"); + Tmax_ = units.convert(T_range[1], "K"); + Pmin_ = units.convert(P_range[0], "Pa"); + Pmax_ = units.convert(P_range[1], "Pa"); + nP_ = coeffs.nColumns(); + nT_ = coeffs.nRows(); + } else { + // ensure that reaction rate can be evaluated (but returns NaN) + coeffs = Array2D(1, 1); + coeffs(0, 0) = NAN; + Tmin_ = 290.; + Tmax_ = 3000.; + Pmin_ = 1.e-7; + Pmax_ = 1.e14; + nP_ = 1; + nT_ = 1; + } + + chebCoeffs_.resize(nP_ * nT_); + dotProd_.resize(nT_); + setup(Tmin_, Tmax_, Pmin_, Pmax_, coeffs); +} + +void Chebyshev::setup(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs) { double logPmin = std::log10(Pmin); double logPmax = std::log10(Pmax); @@ -238,6 +345,41 @@ ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, } } +void Chebyshev::getParameters(AnyMap& rateNode, const Units& rate_units) const +{ + double A = chebCoeffs_[0]; + if (std::isnan(A)) { + // Return empty/unmodified AnyMap + return; + } + rateNode["temperature-range"].setQuantity({Tmin(), Tmax()}, "K"); + rateNode["pressure-range"].setQuantity({Pmin(), Pmax()}, "Pa"); + const auto& coeffs1d = coeffs(); + size_t nT = nTemperature(); + size_t nP = nPressure(); + std::vector coeffs2d(nT, vector_fp(nP)); + for (size_t i = 0; i < nT; i++) { + for (size_t j = 0; j < nP; j++) { + coeffs2d[i][j] = coeffs1d[nP*i + j]; + } + } + // Unit conversions must take place later, after the destination unit system + // is known. A lambda function is used here to override the default behavior + Units rate_units2 = rate_units; + auto converter = [rate_units2](AnyValue& coeffs, const UnitSystem& units) { + if (rate_units2.factor() != 0.0) { + coeffs.asVector()[0][0] += std::log10(units.convertFrom(1.0, rate_units2)); + } else if (units.getDelta(UnitSystem()).size()) { + throw CanteraError("Chebyshev::getParameters lambda", + "Cannot convert rate constant with unknown dimensions to a " + "non-default unit system"); + } + }; + AnyValue coeffs; + coeffs = std::move(coeffs2d); + rateNode["data"].setQuantity(coeffs, converter); +} + BMSurfaceArrhenius::BMSurfaceArrhenius() : m_b(0.0) , m_A(0.0) @@ -272,4 +414,19 @@ void BMSurfaceArrhenius::addCoverageDependence(size_t k, double a, } } +ChebyshevRate::ChebyshevRate() + : Chebyshev() +{ + warn_deprecated("ChebyshevRate::ChebyshevRate", + "Renamed to Chebyshev. Behavior will change after Cantera 2.6. " + "For future behavior, refer to ChebyshevRate3"); +} + +ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs) + : Chebyshev(Tmin, Tmax, Pmin, Pmax, coeffs) +{ + warn_deprecated("ChebyshevRate::ChebyshevRate", + "Renamed to Chebyshev. Behavior will change after Cantera 2.6. " + "For future behavior, refer to ChebyshevRate3");} } diff --git a/test/data/kineticsfromscratch.yaml b/test/data/kineticsfromscratch.yaml index 92cced86289..f59c2c31559 100644 --- a/test/data/kineticsfromscratch.yaml +++ b/test/data/kineticsfromscratch.yaml @@ -45,3 +45,5 @@ reactions: - [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] +- equation: H + 2 O2 <=> HO2 + O2 # Reaction 6 + rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0} diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index 42845753a3d..291af13dda9 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -57,7 +57,7 @@ TEST_F(KineticsFromScratch, add_elementary_reaction) Composition reac = parseCompString("O:1 H2:1"); Composition prod = parseCompString("H:1 OH:1"); Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); kin.addReaction(R); check_rates(0); @@ -73,7 +73,7 @@ TEST_F(KineticsFromScratch, add_three_body_reaction) Arrhenius 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 R = make_shared(reac, prod, rate, tbody); kin.addReaction(R); check_rates(1); @@ -86,7 +86,7 @@ TEST_F(KineticsFromScratch, undefined_third_body) Arrhenius 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 R = make_shared(reac, prod, rate, tbody); ASSERT_THROW(kin.addReaction(R), CanteraError); } @@ -98,7 +98,7 @@ TEST_F(KineticsFromScratch, skip_undefined_third_body) Arrhenius 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 R = make_shared(reac, prod, rate, tbody); kin.skipUndeclaredThirdBodies(true); kin.addReaction(R); @@ -144,7 +144,7 @@ TEST_F(KineticsFromScratch, add_plog_reaction) { 100.0*101325, Arrhenius(5.963200e+56, -11.529, 52599.6 / GasConst_cal_mol_K) } }; - auto R = make_shared(reac, prod, Plog(rates)); + auto R = make_shared(reac, prod, Plog(rates)); kin.addReaction(R); check_rates(3); } @@ -160,7 +160,7 @@ TEST_F(KineticsFromScratch, plog_invalid_rate) { 100.0*101325, Arrhenius(5.9632e+56, -11.529, 52599.6 / GasConst_cal_mol_K) } }; - auto R = make_shared(reac, prod, Plog(rates)); + auto R = make_shared(reac, prod, Plog(rates)); ASSERT_THROW(kin.addReaction(R), CanteraError); } @@ -189,9 +189,9 @@ TEST_F(KineticsFromScratch, add_chebyshev_reaction) coeffs(2,1) = 2.6889e-01; coeffs(2,2) = 9.4806e-02; coeffs(2,3) = -7.6385e-03; - ChebyshevRate rate(290, 3000, 1000.0, 10000000.0, coeffs); + Chebyshev rate(290, 3000, 1000.0, 10000000.0, coeffs); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); kin.addReaction(R); check_rates(4); } @@ -201,7 +201,7 @@ TEST_F(KineticsFromScratch, undeclared_species) Composition reac = parseCompString("CO:1 OH:1"); Composition prod = parseCompString("CO2:1 H:1"); Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); ASSERT_THROW(kin.addReaction(R), CanteraError); ASSERT_EQ((size_t) 0, kin.nReactions()); @@ -212,7 +212,7 @@ TEST_F(KineticsFromScratch, skip_undeclared_species) Composition reac = parseCompString("CO:1 OH:1"); Composition prod = parseCompString("CO2:1 H:1"); Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); kin.skipUndeclaredSpecies(true); kin.addReaction(R); @@ -224,7 +224,7 @@ TEST_F(KineticsFromScratch, negative_A_error) Composition reac = parseCompString("O:1 H2:1"); Composition prod = parseCompString("H:1 OH:1"); Arrhenius rate(-3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); ASSERT_THROW(kin.addReaction(R), CanteraError); ASSERT_EQ((size_t) 0, kin.nReactions()); @@ -235,7 +235,7 @@ TEST_F(KineticsFromScratch, allow_negative_A) Composition reac = parseCompString("O:1 H2:1"); Composition prod = parseCompString("H:1 OH:1"); Arrhenius rate(-3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); R->allow_negative_pre_exponential_factor = true; kin.addReaction(R); @@ -247,7 +247,7 @@ TEST_F(KineticsFromScratch, invalid_reversible_with_orders) Composition reac = parseCompString("O:1 H2:1"); Composition prod = parseCompString("H:1 OH:1"); Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); R->orders["H2"] = 0.5; ASSERT_THROW(kin.addReaction(R), CanteraError); @@ -259,7 +259,7 @@ TEST_F(KineticsFromScratch, negative_order_override) Composition reac = parseCompString("O:1 H2:1"); Composition prod = parseCompString("H:1 OH:1"); Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); R->reversible = false; R->allow_negative_orders = true; R->orders["H2"] = - 0.5; @@ -273,7 +273,7 @@ TEST_F(KineticsFromScratch, invalid_negative_orders) Composition reac = parseCompString("O:1 H2:1"); Composition prod = parseCompString("H:1 OH:1"); Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); R->reversible = false; R->orders["H2"] = - 0.5; @@ -286,7 +286,7 @@ TEST_F(KineticsFromScratch, nonreactant_order_override) Composition reac = parseCompString("O:1 H2:1"); Composition prod = parseCompString("H:1 OH:1"); Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); R->reversible = false; R->allow_nonreactant_orders = true; R->orders["OH"] = 0.5; @@ -300,7 +300,7 @@ TEST_F(KineticsFromScratch, invalid_nonreactant_order) Composition reac = parseCompString("O:1 H2:1"); Composition prod = parseCompString("H:1 OH:1"); Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); - auto R = make_shared(reac, prod, rate); + auto R = make_shared(reac, prod, rate); R->reversible = false; R->orders["OH"] = 0.5; diff --git a/test/kinetics/kineticsFromScratch3.cpp b/test/kinetics/kineticsFromScratch3.cpp new file mode 100644 index 00000000000..603a2d24143 --- /dev/null +++ b/test/kinetics/kineticsFromScratch3.cpp @@ -0,0 +1,407 @@ +#include "gtest/gtest.h" +#include "cantera/kinetics/KineticsFactory.h" +#include "cantera/thermo/IdealGasPhase.h" +#include "cantera/thermo/SurfPhase.h" +#include "cantera/thermo/Species.h" +#include "cantera/kinetics/GasKinetics.h" +#include "cantera/kinetics/InterfaceKinetics.h" +#include "cantera/kinetics/FalloffFactory.h" +#include "cantera/base/Array.h" +#include "cantera/base/stringUtils.h" + +using namespace Cantera; + +class KineticsFromScratch3 : public testing::Test +{ +public: + KineticsFromScratch3() + : p("../data/kineticsfromscratch.yaml") + , p_ref("../data/kineticsfromscratch.yaml") + { + std::vector th; + th.push_back(&p_ref); + kin_ref = newKinetics(th, "../data/kineticsfromscratch.yaml", "ohmech"); + kin.addPhase(p); + kin.init(); + } + + IdealGasPhase p; + IdealGasPhase p_ref; + GasKinetics kin; + unique_ptr kin_ref; + + //! iRef is the index of the corresponding reaction in the reference mech + void check_rates(int iRef) { + ASSERT_EQ((size_t) 1, kin.nReactions()); + + std::string X = "O:0.02 H2:0.2 O2:0.5 H:0.03 OH:0.05 H2O:0.1 HO2:0.01"; + p.setState_TPX(1200, 5*OneAtm, X); + p_ref.setState_TPX(1200, 5*OneAtm, X); + + vector_fp k(1), k_ref(kin_ref->nReactions()); + + kin.getFwdRateConstants(&k[0]); + kin_ref->getFwdRateConstants(&k_ref[0]); + EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + + kin.getRevRateConstants(&k[0]); + kin_ref->getRevRateConstants(&k_ref[0]); + EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + } +}; + +TEST_F(KineticsFromScratch3, add_elementary_reaction) +{ + // reaction 0: + // reaction('O + H2 <=> H + OH', [3.870000e+01, 2.7, 6260.0]) + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + ArrheniusRate rate(3.87e1, 2.7, 2.619184e+07); + auto R = make_shared(reac, prod, rate); + + kin.addReaction(R); + check_rates(0); +} + +TEST_F(KineticsFromScratch3, add_three_body_reaction) +{ + // reaction 1: + // three_body_reaction('2 O + M <=> O2 + M', [1.200000e+11, -1.0, 0.0], + // 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); + + kin.addReaction(R); + check_rates(1); +} + +TEST_F(KineticsFromScratch3, 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); + + ASSERT_THROW(kin.addReaction(R), CanteraError); +} + +TEST_F(KineticsFromScratch3, 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); + + kin.skipUndeclaredThirdBodies(true); + kin.addReaction(R); + ASSERT_EQ((size_t) 1, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, add_plog_reaction) +{ + // reaction 3: + // pdep_arrhenius('H2 + O2 <=> 2 OH', + // [(0.01, 'atm'), 1.212400e+16, -0.5779, 10872.7], + // [(1.0, 'atm'), 4.910800e+31, -4.8507, 24772.8], + // [(10.0, 'atm'), 1.286600e+47, -9.0246, 39796.5], + // [(100.0, 'atm'), 5.963200e+56, -11.529, 52599.6]) + Composition reac = parseCompString("H2:1, O2:1"); + Composition prod = parseCompString("OH:2"); + std::multimap rates { + { 0.01*101325, Arrhenius(1.212400e+16, -0.5779, 10872.7 / GasConst_cal_mol_K) }, + { 1.0*101325, Arrhenius(4.910800e+31, -4.8507, 24772.8 / GasConst_cal_mol_K) }, + { 10.0*101325, Arrhenius(1.286600e+47, -9.0246, 39796.5 / GasConst_cal_mol_K) }, + { 100.0*101325, Arrhenius(5.963200e+56, -11.529, 52599.6 / GasConst_cal_mol_K) } + }; + + auto R = make_shared(reac, prod, PlogRate(rates)); + kin.addReaction(R); + check_rates(3); +} + +TEST_F(KineticsFromScratch3, plog_invalid_rate) +{ + Composition reac = parseCompString("H2:1, O2:1"); + Composition prod = parseCompString("OH:2"); + std::multimap rates { + { 0.01*101325, Arrhenius(1.2124e+16, -0.5779, 10872.7 / GasConst_cal_mol_K) }, + { 10.0*101325, Arrhenius(1e15, -1, 10000 / GasConst_cal_mol_K) }, + { 10.0*101325, Arrhenius(-2e20, -2.0, 20000 / GasConst_cal_mol_K) }, + { 100.0*101325, Arrhenius(5.9632e+56, -11.529, 52599.6 / GasConst_cal_mol_K) } + }; + + auto R = make_shared(reac, prod, PlogRate(rates)); + ASSERT_THROW(kin.addReaction(R), CanteraError); +} + +TEST_F(KineticsFromScratch3, add_chebyshev_reaction) +{ + // reaction 4: + // chebyshev_reaction( + // 'HO2 <=> OH + O', + // Tmin=290.0, Tmax=3000.0, + // Pmin=(0.0098692326671601278, 'atm'), Pmax=(98.692326671601279, 'atm'), + // coeffs=[[ 8.2883e+00, -1.1397e+00, -1.2059e-01, 1.6034e-02], + // [ 1.9764e+00, 1.0037e+00, 7.2865e-03, -3.0432e-02], + // [ 3.1770e-01, 2.6889e-01, 9.4806e-02, -7.6385e-03]]) + Composition reac = parseCompString("HO2:1"); + Composition prod = parseCompString("OH:1 O:1"); + Array2D coeffs(3, 4); + coeffs(0,0) = 8.2883e+00; + coeffs(0,1) = -1.1397e+00; + coeffs(0,2) = -1.2059e-01; + coeffs(0,3) = 1.6034e-02; + coeffs(1,0) = 1.9764e+00; + coeffs(1,1) = 1.0037e+00; + coeffs(1,2) = 7.2865e-03; + coeffs(1,3) = -3.0432e-02; + coeffs(2,0) = 3.1770e-01; + coeffs(2,1) = 2.6889e-01; + coeffs(2,2) = 9.4806e-02; + coeffs(2,3) = -7.6385e-03; + ChebyshevRate3 rate(290, 3000, 1000.0, 10000000.0, coeffs); + + auto R = make_shared(reac, prod, rate); + kin.addReaction(R); + check_rates(4); +} + +TEST_F(KineticsFromScratch3, undeclared_species) +{ + Composition reac = parseCompString("CO:1 OH:1"); + Composition prod = parseCompString("CO2:1 H:1"); + ArrheniusRate rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + + ASSERT_THROW(kin.addReaction(R), CanteraError); + ASSERT_EQ((size_t) 0, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, skip_undeclared_species) +{ + Composition reac = parseCompString("CO:1 OH:1"); + Composition prod = parseCompString("CO2:1 H:1"); + ArrheniusRate rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + + kin.skipUndeclaredSpecies(true); + kin.addReaction(R); + ASSERT_EQ((size_t) 0, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, negative_A_error) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + ArrheniusRate rate(-3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + + ASSERT_THROW(kin.addReaction(R), CanteraError); + ASSERT_EQ((size_t) 0, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, allow_negative_A) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + ArrheniusRate rate(-3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + auto rr = std::dynamic_pointer_cast(R->rate()); + rr->allow_negative_pre_exponential_factor = true; + + kin.addReaction(R); + ASSERT_EQ((size_t) 1, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, invalid_reversible_with_orders) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + ArrheniusRate rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + R->orders["H2"] = 0.5; + + ASSERT_THROW(kin.addReaction(R), CanteraError); + ASSERT_EQ((size_t) 0, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, negative_order_override) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + ArrheniusRate rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + R->reversible = false; + R->allow_negative_orders = true; + R->orders["H2"] = - 0.5; + + kin.addReaction(R); + ASSERT_EQ((size_t) 1, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, invalid_negative_orders) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + ArrheniusRate rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + R->reversible = false; + R->orders["H2"] = - 0.5; + + ASSERT_THROW(kin.addReaction(R), CanteraError); + ASSERT_EQ((size_t) 0, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, nonreactant_order_override) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + ArrheniusRate rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + R->reversible = false; + R->allow_nonreactant_orders = true; + R->orders["OH"] = 0.5; + + kin.addReaction(R); + ASSERT_EQ((size_t) 1, kin.nReactions()); +} + +TEST_F(KineticsFromScratch3, invalid_nonreactant_order) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + ArrheniusRate rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + auto R = make_shared(reac, prod, rate); + R->reversible = false; + R->orders["OH"] = 0.5; + + ASSERT_THROW(kin.addReaction(R), CanteraError); + ASSERT_EQ((size_t) 0, kin.nReactions()); +} + +class KineticsAddSpecies3 : public testing::Test +{ +public: + KineticsAddSpecies3() + : p_ref("../data/kineticsfromscratch.yaml") + { + std::vector th; + th.push_back(&p_ref); + kin_ref = newKinetics(th, "../data/kineticsfromscratch.yaml", "ohmech"); + kin.addPhase(p); + + std::vector> S = getSpecies( + AnyMap::fromYamlFile("h2o2.yaml")["species"]); + for (auto sp : S) { + species[sp->name] = sp; + } + reactions = getReactions( + AnyMap::fromYamlFile("../data/kineticsfromscratch.yaml")["reactions"], + *kin_ref); + } + + IdealGasPhase p; + IdealGasPhase p_ref; + GasKinetics kin; + unique_ptr kin_ref; + std::vector> reactions; + std::map> species; + + void check_rates(size_t N, const std::string& X) { + for (size_t i = 0; i < kin_ref->nReactions(); i++) { + if (i >= N) { + kin_ref->setMultiplier(i, 0); + } else { + kin_ref->setMultiplier(i, 1); + } + } + p.setState_TPX(1200, 5*OneAtm, X); + p_ref.setState_TPX(1200, 5*OneAtm, X); + + vector_fp k(kin.nReactions()), k_ref(kin_ref->nReactions()); + vector_fp w(kin.nTotalSpecies()), w_ref(kin_ref->nTotalSpecies()); + + kin.getFwdRateConstants(k.data()); + kin_ref->getFwdRateConstants(k_ref.data()); + for (size_t i = 0; i < kin.nReactions(); i++) { + EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + } + + kin.getFwdRatesOfProgress(k.data()); + kin_ref->getFwdRatesOfProgress(k_ref.data()); + for (size_t i = 0; i < kin.nReactions(); i++) { + EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + } + + kin.getRevRateConstants(k.data()); + kin_ref->getRevRateConstants(k_ref.data()); + for (size_t i = 0; i < kin.nReactions(); i++) { + EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + } + + kin.getRevRatesOfProgress(k.data()); + kin_ref->getRevRatesOfProgress(k_ref.data()); + for (size_t i = 0; i < kin.nReactions(); i++) { + EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + } + + kin.getCreationRates(w.data()); + kin_ref->getCreationRates(w_ref.data()); + for (size_t i = 0; i < kin.nTotalSpecies(); i++) { + size_t iref = p_ref.speciesIndex(p.speciesName(i)); + EXPECT_DOUBLE_EQ(w_ref[iref], w[i]) << "sp = " << p.speciesName(i) << "; N = " << N; + } + } +}; + +TEST_F(KineticsAddSpecies3, add_species_sequential) +{ + ASSERT_EQ((size_t) 0, kin.nReactions()); + + for (auto s : {"AR", "O", "H2", "H", "OH"}) { + p.addSpecies(species[s]); + } + kin.addReaction(reactions[0]); + ASSERT_EQ(5, (int) kin.nTotalSpecies()); + check_rates(1, "O:0.001, H2:0.1, H:0.005, OH:0.02, AR:0.88"); + + p.addSpecies(species["O2"]); + p.addSpecies(species["H2O"]); + kin.addReaction(reactions[1]); + ASSERT_EQ(7, (int) kin.nTotalSpecies()); + ASSERT_EQ(2, (int) kin.nReactions()); + check_rates(2, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.38"); + + p.addSpecies(species["H2O2"]); + kin.addReaction(reactions[2]); + kin.addReaction(reactions[3]); + check_rates(4, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.38"); // no change + check_rates(4, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.35, H2O2:0.03"); + + p.addSpecies(species["HO2"]); + kin.addReaction(reactions[4]); + check_rates(5, "O:0.01, H2:0.1, H:0.02, OH:0.03, O2:0.4, AR:0.3, H2O2:0.03, HO2:0.01"); +} + +TEST_F(KineticsAddSpecies3, add_species_err_first) +{ + for (auto s : {"AR", "O", "H2", "H"}) { + p.addSpecies(species[s]); + } + ASSERT_THROW(kin.addReaction(reactions[0]), CanteraError); + ASSERT_EQ((size_t) 0, kin.nReactions()); + + p.addSpecies(species["OH"]); + kin.addReaction(reactions[0]); + ASSERT_EQ(5, (int) kin.nTotalSpecies()); + ASSERT_EQ((size_t) 1, kin.nReactions()); + check_rates(1, "O:0.001, H2:0.1, H:0.005, OH:0.02, AR:0.88"); +} diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 261909768da..93fb086c53c 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -10,7 +10,24 @@ using namespace Cantera; -TEST(Reaction, ElementaryFromYaml) +TEST(ReactionRate, ModifyArrheniusRate) +{ + auto sol = newSolution("gri30.yaml", "", "None"); + AnyMap rxn = AnyMap::fromYamlString( + "{equation: N + NO <=> N2 + O," + " rate-constant: [-2.70000E+13 cm^3/mol/s, 0, 355 cal/mol]," + " negative-A: true}"); + + auto R = newReaction(rxn, *(sol->kinetics())); + auto ER = dynamic_cast(*R); + + auto rr = std::dynamic_pointer_cast(ER.rate()); + EXPECT_TRUE(rr->allow_negative_pre_exponential_factor); + rr->allow_negative_pre_exponential_factor = false; + EXPECT_FALSE(rr->allow_negative_pre_exponential_factor); +} + +TEST(Reaction, ElementaryFromYaml3) { auto sol = newSolution("gri30.yaml", "", "None"); AnyMap rxn = AnyMap::fromYamlString( @@ -22,15 +39,36 @@ TEST(Reaction, ElementaryFromYaml) EXPECT_EQ(R->reactants.at("NO"), 1); EXPECT_EQ(R->products.at("N2"), 1); EXPECT_EQ(R->type(), "elementary"); + EXPECT_FALSE(R->allow_negative_orders); + + const auto& rate = std::dynamic_pointer_cast(R->rate()); + EXPECT_TRUE(rate->allow_negative_pre_exponential_factor); + EXPECT_DOUBLE_EQ(rate->preExponentialFactor(), -2.7e10); + EXPECT_DOUBLE_EQ(rate->activationEnergy_R(), 355 / GasConst_cal_mol_K); +} + +TEST(Reaction, ElementaryFromYaml2) +{ + auto sol = newSolution("gri30.yaml"); + AnyMap rxn = AnyMap::fromYamlString( + "{equation: N + NO <=> N2 + O," + " type: elementary-legacy," + " rate-constant: [-2.70000E+13 cm^3/mol/s, 0, 355 cal/mol]," + " negative-A: true}"); + + auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_EQ(R->reactants.at("NO"), 1); + EXPECT_EQ(R->products.at("N2"), 1); + EXPECT_EQ(R->type(), "elementary-legacy"); - auto ER = dynamic_cast(*R); + auto ER = dynamic_cast(*R); EXPECT_DOUBLE_EQ(ER.rate.preExponentialFactor(), -2.7e10); EXPECT_DOUBLE_EQ(ER.rate.activationEnergy_R(), 355 / GasConst_cal_mol_K); EXPECT_TRUE(ER.allow_negative_pre_exponential_factor); EXPECT_FALSE(ER.allow_negative_orders); } -TEST(Reaction, ThreeBodyFromYaml1) +TEST(Reaction, ThreeBodyFromYaml3) { auto sol = newSolution("gri30.yaml", "", "None"); AnyMap rxn = AnyMap::fromYamlString( @@ -42,14 +80,14 @@ TEST(Reaction, ThreeBodyFromYaml1) auto R = newReaction(rxn, *(sol->kinetics())); EXPECT_EQ(R->reactants.count("M"), (size_t) 0); EXPECT_EQ(R->type(), "three-body"); + EXPECT_DOUBLE_EQ(R->thirdBody()->efficiencies["H2O"], 5.0); + EXPECT_DOUBLE_EQ(R->thirdBody()->default_efficiency, 1.0); - auto TBR = dynamic_cast(*R); - EXPECT_DOUBLE_EQ(TBR.rate.preExponentialFactor(), 1.2e11); - EXPECT_DOUBLE_EQ(TBR.third_body.efficiencies["H2O"], 5.0); - EXPECT_DOUBLE_EQ(TBR.third_body.default_efficiency, 1.0); + const auto& rate = std::dynamic_pointer_cast(R->rate()); + EXPECT_DOUBLE_EQ(rate->preExponentialFactor(), 1.2e11); } -TEST(Reaction, ThreeBodyFromYaml2) +TEST(Reaction, ThreeBodyFromYamlMissingM) { auto sol = newSolution("gri30.yaml", "", "None"); AnyMap rxn = AnyMap::fromYamlString( @@ -60,6 +98,25 @@ TEST(Reaction, ThreeBodyFromYaml2) EXPECT_THROW(newReaction(rxn, *(sol->kinetics())), CanteraError); } +TEST(Reaction, ThreeBodyFromYaml2) +{ + auto sol = newSolution("gri30.yaml"); + AnyMap rxn = AnyMap::fromYamlString( + "{equation: 2 O + M = O2 + M," + " type: three-body-legacy," + " 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_EQ(R->reactants.count("M"), (size_t) 0); + EXPECT_EQ(R->type(), "three-body-legacy"); + + auto TBR = dynamic_cast(*R); + EXPECT_DOUBLE_EQ(TBR.rate.preExponentialFactor(), 1.2e11); + EXPECT_DOUBLE_EQ(TBR.third_body.efficiencies["H2O"], 5.0); + EXPECT_DOUBLE_EQ(TBR.third_body.default_efficiency, 1.0); +} + TEST(Reaction, FalloffFromYaml1) { auto sol = newSolution("gri30.yaml", "", "None"); @@ -134,8 +191,7 @@ TEST(Reaction, PlogFromYaml) "- {P: 1.01325 MPa, A: 1.680000e+16, b: -0.6, Ea: 14754.0}"); auto R = newReaction(rxn, *(sol->kinetics())); - auto PR = dynamic_cast(*R); - const auto& rates = PR.rate.rates(); + const auto& rates = std::dynamic_pointer_cast(R->rate())->rates(); EXPECT_EQ(rates.size(), (size_t) 4); EXPECT_NEAR(rates[0].first, 0.039474 * OneAtm, 1e-6); EXPECT_NEAR(rates[2].first, OneAtm, 1e-6); @@ -162,15 +218,15 @@ TEST(Reaction, ChebyshevFromYaml) auto R = newReaction(rxn, *(sol->kinetics())); EXPECT_EQ(R->reactants.size(), (size_t) 1); - auto CR = dynamic_cast(*R); + const auto rate = std::dynamic_pointer_cast(R->rate()); double logP = std::log10(2e6); double T = 1800; - CR.rate.update_C(&logP); - EXPECT_EQ(CR.rate.nTemperature(), (size_t) 6); - EXPECT_EQ(CR.rate.nPressure(), (size_t) 4); - EXPECT_DOUBLE_EQ(CR.rate.Tmax(), 3000); - EXPECT_DOUBLE_EQ(CR.rate.Pmin(), 1000); - EXPECT_NEAR(CR.rate.updateRC(std::log(T), 1.0/T), 130512.2773948636, 1e-9); + rate->update_C(&logP); + EXPECT_EQ(rate->nTemperature(), (size_t) 6); + EXPECT_EQ(rate->nPressure(), (size_t) 4); + EXPECT_DOUBLE_EQ(rate->Tmax(), 3000); + EXPECT_DOUBLE_EQ(rate->Pmin(), 1000); + EXPECT_NEAR(rate->updateRC(std::log(T), 1.0/T), 130512.2773948636, 1e-9); } TEST(Reaction, BlowersMaselFromYaml) @@ -218,8 +274,9 @@ TEST(Kinetics, GasKineticsFromYaml1) EXPECT_EQ(R->reactants.at("NO"), 1); EXPECT_EQ(R->products.at("N2"), 1); EXPECT_EQ(R->id, "NOx-R1"); - const auto& ER = std::dynamic_pointer_cast(R); - EXPECT_DOUBLE_EQ(ER->rate.preExponentialFactor(), 2.7e10); + const auto& ER = std::dynamic_pointer_cast(R); + const auto& rate = std::dynamic_pointer_cast(ER->rate()); + EXPECT_DOUBLE_EQ(rate->preExponentialFactor(), 2.7e10); } TEST(Kinetics, GasKineticsFromYaml2) @@ -400,7 +457,6 @@ class ReactionToYaml : public testing::Test void compareReactions() { auto kin = soln->kinetics(); EXPECT_EQ(kin->reactionString(iOld), kin->reactionString(iNew)); - EXPECT_EQ(kin->reactionTypeStr(iOld), kin->reactionTypeStr(iNew)); EXPECT_EQ(kin->isReversible(iOld), kin->isReversible(iNew)); vector_fp kf(kin->nReactions()), kr(kin->nReactions()); @@ -427,7 +483,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_TRUE(std::dynamic_pointer_cast(duplicate)); + EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); compareReactions(); } @@ -436,7 +492,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)); compareReactions(); } @@ -474,7 +530,7 @@ TEST_F(ReactionToYaml, pdepArrhenius) soln = newSolution("pdep-test.yaml"); soln->thermo()->setState_TPY(1000, 2e5, "R2:1, H:0.1, P2A:2, P2B:0.3"); duplicateReaction(1); - EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); compareReactions(); soln->thermo()->setState_TPY(1100, 1e3, "R2:1, H:0.2, P2A:2, P2B:0.3"); compareReactions(); @@ -485,7 +541,7 @@ TEST_F(ReactionToYaml, Chebyshev) soln = newSolution("pdep-test.yaml"); soln->thermo()->setState_TPY(1000, 2e5, "R6:1, P6A:2, P6B:0.3"); duplicateReaction(5); - EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); compareReactions(); } @@ -524,9 +580,9 @@ TEST_F(ReactionToYaml, electrochemical) TEST_F(ReactionToYaml, unconvertible1) { - ElementaryReaction R({{"H2", 1}, {"OH", 1}}, - {{"H2O", 1}, {"H", 1}}, - Arrhenius(1e5, -1.0, 12.5)); + ElementaryReaction2 R({{"H2", 1}, {"OH", 1}}, + {{"H2O", 1}, {"H", 1}}, + Arrhenius(1e5, -1.0, 12.5)); AnyMap params = R.parameters(); UnitSystem U{"g", "cm", "mol"}; params.setUnits(U); @@ -536,9 +592,9 @@ TEST_F(ReactionToYaml, unconvertible1) TEST_F(ReactionToYaml, unconvertible2) { Array2D coeffs(2, 2, 1.0); - ChebyshevReaction R({{"H2", 1}, {"OH", 1}}, - {{"H2O", 1}, {"H", 1}}, - ChebyshevRate(273, 3000, 1e2, 1e7, coeffs)); + ChebyshevReaction2 R({{"H2", 1}, {"OH", 1}}, + {{"H2O", 1}, {"H", 1}}, + Chebyshev(273, 3000, 1e2, 1e7, coeffs)); UnitSystem U{"g", "cm", "mol"}; AnyMap params = R.parameters(); params.setUnits(U);