diff --git a/SConstruct b/SConstruct index da7b332a958..25ad07f4698 100644 --- a/SConstruct +++ b/SConstruct @@ -1746,6 +1746,7 @@ def postBuildMessage(target, source, env): print("*******************************************************") print("Compilation completed successfully.\n") print("- To run the test suite, type 'scons test'.") + print("- To list available tests, type 'scons test-help'.") if env['googletest'] == 'none': print(" WARNING: You set the 'googletest' to 'none' and all it's tests will be skipped.") if os.name == 'nt': diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h index 52672639c09..a778e199e15 100644 --- a/include/cantera/kinetics/BulkKinetics.h +++ b/include/cantera/kinetics/BulkKinetics.h @@ -11,6 +11,7 @@ #include "Kinetics.h" #include "RateCoeffMgr.h" +#include "cantera/kinetics/MultiRate.h" namespace Cantera { @@ -37,6 +38,8 @@ class BulkKinetics : public Kinetics bool doIrreversible = false); virtual bool addReaction(shared_ptr r); + virtual void modifyReaction(size_t i, shared_ptr rNew); + virtual void resizeSpecies(); virtual void setMultiplier(size_t i, double f); @@ -46,6 +49,10 @@ class BulkKinetics : public Kinetics virtual void addElementaryReaction(ElementaryReaction& r); virtual void modifyElementaryReaction(size_t i, ElementaryReaction& rNew); + //! Vector of rate handlers + std::vector> m_bulk_rates; + std::map m_bulk_types; //!< Mapping of rate handlers + Rate1 m_rates; std::vector m_revindex; //!< Indices of reversible reactions std::vector m_irrev; //!< Indices of irreversible reactions diff --git a/include/cantera/kinetics/FalloffMgr.h b/include/cantera/kinetics/FalloffMgr.h index f665cabe89c..a89f97ff049 100644 --- a/include/cantera/kinetics/FalloffMgr.h +++ b/include/cantera/kinetics/FalloffMgr.h @@ -9,6 +9,7 @@ #define CT_FALLOFFMGR_H #include "Falloff.h" +#include "cantera/base/global.h" namespace Cantera { @@ -32,13 +33,35 @@ class FalloffMgr * determine which array entry is modified in method pr_to_falloff. * @param reactionType Either `FALLOFF_RXN` or `CHEMACT_RXN` * @param f The falloff function. + * + * @deprecated To be removed after Cantera 2.6. */ void install(size_t rxn, int reactionType, shared_ptr f) { + warn_deprecated("FalloffMgr::install()", + "To be removed after Cantera 2.6. Specify reaction type using " + "string instead."); + + m_rxn.push_back(rxn); + m_offset.push_back(m_worksize); + m_worksize += f->workSize(); + m_falloff.push_back(f); + m_isfalloff.push_back(reactionType == FALLOFF_RXN); + m_indices[rxn] = m_falloff.size()-1; + } + + //! Install a new falloff function calculator. + /* + * @param rxn Index of the falloff reaction. This will be used to + * determine which array entry is modified in method pr_to_falloff. + * @param type Reaction type identifier. + * @param f The falloff function. + */ + void install(size_t rxn, std::string type, shared_ptr f) { m_rxn.push_back(rxn); m_offset.push_back(m_worksize); m_worksize += f->workSize(); m_falloff.push_back(f); - m_reactionType.push_back(reactionType); + m_isfalloff.push_back(type == "falloff"); m_indices[rxn] = m_falloff.size()-1; } @@ -76,14 +99,14 @@ class FalloffMgr void pr_to_falloff(doublereal* values, const doublereal* work) { for (size_t i = 0; i < m_rxn.size(); i++) { double pr = values[m_rxn[i]]; - if (m_reactionType[i] == FALLOFF_RXN) { + if (m_isfalloff[i]) { // Pr / (1 + Pr) * F values[m_rxn[i]] *= - m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr); + m_falloff[i]->F(pr, work + m_offset[i]) / (1.0 + pr); } else { // 1 / (1 + Pr) * F values[m_rxn[i]] = - m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr); + m_falloff[i]->F(pr, work + m_offset[i]) / (1.0 + pr); } } } @@ -96,7 +119,7 @@ class FalloffMgr size_t m_worksize; //! Distinguish between falloff and chemically activated reactions - vector_int m_reactionType; + std::vector m_isfalloff; //! map of external reaction index to local index std::map m_indices; diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index ff2c1ee7554..f448243f937 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -13,6 +13,7 @@ #include "StoichManager.h" #include "cantera/base/ValueCache.h" +#include "cantera/kinetics/ReactionFactory.h" namespace Cantera { @@ -601,9 +602,18 @@ class Kinetics * are specific to the particular kinetics manager. * * @param i reaction index + * + * @deprecated To be changed after Cantera 2.6. */ virtual int reactionType(size_t i) const; + /** + * String specifying the type of reaction. + * + * @param i reaction index + */ + virtual std::string reactionTypeStr(size_t i) const; + /** * True if reaction i has been declared to be reversible. If isReversible(i) * is false, then the reverse rate of progress for reaction i is always diff --git a/include/cantera/kinetics/MultiRate.h b/include/cantera/kinetics/MultiRate.h new file mode 100644 index 00000000000..e2fad2002d1 --- /dev/null +++ b/include/cantera/kinetics/MultiRate.h @@ -0,0 +1,127 @@ +/** + * @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 +// at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_MULTIRATE_H +#define CT_MULTIRATE_H + +#include "cantera/kinetics/ReactionRate.h" + +namespace Cantera +{ + +//! An abstract base class for evaluating all reactions of a particular type. +/** + * Because this class has no template parameters, the `Kinetics` object + * can store all of these rate coefficient evaluators as a + * `vector>`. + * + * @todo At the moment, implemented methods are specific to `BulkKinetics`, + * 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 +{ +public: + virtual ~MultiRateBase() {} + + //! Add reaction rate object to the evaluator + //! @param rxn_index index of reaction + //! @param rate reaction rate object + virtual void add(const size_t rxn_index, + ReactionRateBase& rate) = 0; + + //! Replace reaction rate object handled by the evaluator + //! @param rxn_index index of reaction + //! @param rate reaction rate object + virtual bool replace(const size_t rxn_index, + ReactionRateBase& rate) = 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; + + //! Update data common to reaction rates of a specific type + //! @param bulk object representing bulk phase + virtual void update(const ThermoPhase& bulk) = 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); + } + + virtual bool replace(const size_t rxn_index, + ReactionRateBase& rate) override { + if (!m_rates.size()) { + throw CanteraError("MultiBulkRate::replace", + "Invalid operation: cannot replace rate object " + "in empty rate handler."); + } else if (typeid(rate) != typeid(RateType)) { + throw CanteraError("MultiBulkRate::replace", + "Invalid operation: cannot replace rate object of type '{}' " + "with a new rate of type '{}'.", + m_rates[0].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); + return true; + } + return false; + } + + 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); + } + } + + virtual void update(const ThermoPhase& bulk) override { + // update common data once for each reaction type + m_shared.update(bulk); + if (RateType::uses_update()) { + // 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); + } + } + } + +protected: + std::vector m_rates; //! Reaction rate objects + std::vector m_rxn; //! Index within overall rate vector + std::map m_indices; //! Mapping of indices + DataType m_shared; +}; + +} + +#endif diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index e8a300f175c..f19b3cf69d7 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -9,7 +9,7 @@ #define CT_REACTION_H #include "cantera/base/AnyMap.h" -#include "cantera/kinetics/RxnRates.h" +#include "cantera/kinetics/ReactionRate.h" namespace Cantera { @@ -18,12 +18,18 @@ class Kinetics; class Falloff; class XML_Node; -//! Intermediate class which stores data about a reaction and its rate +//! Abstract base class which stores data about a reaction and its rate //! parameterization so that it can be added to a Kinetics object. class Reaction { public: + Reaction(); + Reaction(const Composition& reactants, + const Composition& products); + + //! @deprecated To be removed after Cantera 2.6. explicit Reaction(int type); + //! @deprecated To be removed after Cantera 2.6. Reaction(int type, const Composition& reactants, const Composition& products); virtual ~Reaction() {} @@ -37,12 +43,29 @@ class Reaction //! The chemical equation for this reaction std::string equation() const; + //! The type of reaction + virtual std::string type() const = 0; // pure virtual function + //! Ensure that the rate constant and other parameters for this reaction are //! valid. virtual void validate(); + //! Get validity flag of reaction + bool valid() const { + return m_valid; + } + + //! Set validity flag of reaction + void setValid(bool valid) { + m_valid = valid; + } + //! Type of the reaction. The valid types are listed in the file, //! reaction_defs.h, with constants ending in `RXN`. + /*! + * @deprecated To be removed in Cantera 2.6. + * Superseded by Reaction::type(). + */ int reaction_type; //! Reactant species and stoichiometric coefficients @@ -75,6 +98,35 @@ class Reaction //! Input data used for specific models AnyMap input; + +protected: + //! Flag indicating whether reaction is set up correctly + bool m_valid; +}; + + +//! 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; }; @@ -88,6 +140,10 @@ class ElementaryReaction : public Reaction const Arrhenius& rate); virtual void validate(); + virtual std::string type() const { + return "elementary"; + } + Arrhenius rate; bool allow_negative_pre_exponential_factor; }; @@ -117,6 +173,11 @@ class ThreeBodyReaction : public ElementaryReaction ThreeBodyReaction(); ThreeBodyReaction(const Composition& reactants, const Composition& products, const Arrhenius& rate, const ThirdBody& tbody); + + virtual std::string type() const { + return "three-body"; + } + virtual std::string reactantString() const; virtual std::string productString() const; @@ -134,6 +195,11 @@ class FalloffReaction : public Reaction FalloffReaction(const Composition& reactants, const Composition& products, const Arrhenius& low_rate, const Arrhenius& high_rate, const ThirdBody& tbody); + + virtual std::string type() const { + return "falloff"; + } + virtual std::string reactantString() const; virtual std::string productString() const; virtual void validate(); @@ -165,6 +231,10 @@ class ChemicallyActivatedReaction : public FalloffReaction ChemicallyActivatedReaction(const Composition& reactants, const Composition& products, const Arrhenius& low_rate, const Arrhenius& high_rate, const ThirdBody& tbody); + + virtual std::string type() const { + return "chemically-activated"; + } }; //! A pressure-dependent reaction parameterized by logarithmically interpolating @@ -175,6 +245,11 @@ class PlogReaction : public Reaction PlogReaction(); PlogReaction(const Composition& reactants, const Composition& products, const Plog& rate); + + virtual std::string type() const { + return "pressure-dependent-Arrhenius"; + } + virtual void validate(); Plog rate; }; @@ -188,9 +263,55 @@ class ChebyshevReaction : public Reaction ChebyshevReaction(const Composition& reactants, const Composition& products, const ChebyshevRate& rate); + virtual std::string type() const { + return "Chebyshev"; + } + 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; +}; + + //! Modifications to an InterfaceReaction rate based on a surface species //! coverage. struct CoverageDependency @@ -214,6 +335,10 @@ class InterfaceReaction : public ElementaryReaction InterfaceReaction(const Composition& reactants, const Composition& products, const Arrhenius& rate, bool isStick=false); + virtual std::string type() const { + return "interface"; + } + //! Adjustments to the Arrhenius rate expression dependent on surface //! species coverages. Three coverage parameters (a, E, m) are used for each //! species on which the rate depends. See SurfaceArrhenius for details on @@ -249,15 +374,6 @@ class ElectrochemicalReaction : public InterfaceReaction bool exchange_current_density_formulation; }; -//! Create a new Reaction object for the reaction defined in `rxn_node` -//! -//! @deprecated The XML input format is deprecated and will be removed in -//! Cantera 3.0. -shared_ptr newReaction(const XML_Node& rxn_node); - -//! Create a new Reaction object using the specified parameters -unique_ptr newReaction(const AnyMap& rxn_node, const Kinetics& kin); - //! Create Reaction objects for all `` nodes in an XML document. //! //! The `` nodes are assumed to be children of the `` @@ -286,6 +402,59 @@ std::vector > getReactions(const XML_Node& node); //! Kinetics object `kinetics`. std::vector> getReactions(const AnyValue& items, Kinetics& kinetics); + +//! Parse reaction equation +void parseReactionEquation(Reaction& R, const AnyValue& equation, + const Kinetics& kin); + +//! 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. +//! +//! @todo Rate units will become available as `rate_units` after serialization +//! is implemented. +Units rateCoeffUnits(const Reaction& R, const Kinetics& kin, + int pressure_dependence=0); + +// declarations of setup functions +void setupElementaryReaction(ElementaryReaction&, const XML_Node&); +//! @internal May be changed without notice in future versions +void setupElementaryReaction(ElementaryReaction&, const AnyMap&, + const Kinetics&); + +void setupThreeBodyReaction(ThreeBodyReaction&, const XML_Node&); +//! @internal May be changed without notice in future versions +void setupThreeBodyReaction(ThreeBodyReaction&, const AnyMap&, + const Kinetics&); + +void setupFalloffReaction(FalloffReaction&, const XML_Node&); +//! @internal May be changed without notice in future versions +void setupFalloffReaction(FalloffReaction&, const AnyMap&, + const Kinetics&); + +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 setupChebyshevReaction(ChebyshevReaction&, const XML_Node&); +//! @internal May be changed without notice in future versions +void setupChebyshevReaction(ChebyshevReaction&, const AnyMap&, + const Kinetics&); + +void setupInterfaceReaction(InterfaceReaction&, const XML_Node&); +//! @internal May be changed without notice in future versions +void setupInterfaceReaction(InterfaceReaction&, const AnyMap&, + const Kinetics&); + +void setupElectrochemicalReaction(ElectrochemicalReaction&, + const XML_Node&); +//! @internal May be changed without notice in future versions +void setupElectrochemicalReaction(ElectrochemicalReaction&, + const AnyMap&, const Kinetics&); + } #endif diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h new file mode 100644 index 00000000000..89da9c51e3d --- /dev/null +++ b/include/cantera/kinetics/ReactionData.h @@ -0,0 +1,82 @@ +/** + * @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 +// at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_REACTIONDATA_H +#define CT_REACTIONDATA_H + +namespace Cantera +{ + +class ThermoPhase; + + +//! Data container holding shared data specific to ArrheniusRate +/** + * 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); } + + //! 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); } + + void update(double T) { + m_temperature = T; + m_logT = std::log(T); + m_recipT = 1./T; + } + + void update(const ThermoPhase& bulk); + + double m_temperature; //!< temperature + double m_logT; //!< logarithm of temperature + double m_recipT; //!< inverse of temperature +}; + + +//! Data container holding shared data specific to CustomFunc1Rate +/** + * @warning This class is an experimental part of the %Cantera API and + * may be changed or removed without notice. + */ +struct CustomFunc1Data +{ + CustomFunc1Data() : m_temperature(1.) {} + + //! Constructor based on temperature *T* + CustomFunc1Data(double T) { update(T); } + + //! Constructor based on temperature *T* and pressure *P* + CustomFunc1Data(double T, double P) { update(T); }; + + //! Constructor accessing *bulk* phase definitions + CustomFunc1Data(const ThermoPhase& bulk) { update(bulk); } + + void update(double T) { m_temperature = T; } + + void update(const ThermoPhase& bulk); + + double m_temperature; //!< temperature +}; + +} + +#endif diff --git a/include/cantera/kinetics/ReactionFactory.h b/include/cantera/kinetics/ReactionFactory.h new file mode 100644 index 00000000000..6cf0ff7e09c --- /dev/null +++ b/include/cantera/kinetics/ReactionFactory.h @@ -0,0 +1,144 @@ +/** + * @file ReactionFactory.h + * Factory class for reaction functions. Used by classes + * that implement kinetics + * (see \ref reactionGroup and class \link Cantera::Reaction Reaction\endlink). + */ + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_NEWREACTION_H +#define CT_NEWREACTION_H + +#include "cantera/base/FactoryBase.h" +#include "cantera/kinetics/Reaction.h" + +namespace Cantera +{ + +/** + * Factory class to construct reaction function calculators. + * The reaction factory is accessed through static method factory: + * + * @code + * Reaction* f = ReactionFactory::factory()->newReaction(type, c) + * @endcode + * + * @ingroup reactionGroup + * + * @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 +{ +public: + /** + * Return a pointer to the factory. On the first call, a new instance is + * created. Since there is no need to instantiate more than one factory, + * on all subsequent calls, a pointer to the existing factory is returned. + */ + static ReactionFactory* factory() { + std::unique_lock lock(reaction_mutex); + if (!s_factory) { + s_factory = new ReactionFactory; + } + return s_factory; + } + + virtual void deleteFactory() { + std::unique_lock lock(reaction_mutex); + delete s_factory; + s_factory = 0; + } + + //! 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); + } + } + + //! 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; + } + + //! Set up reaction based on AnyMap node information + /** + * @todo call setup directly from constructor after removal of XML support. + */ + 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); + } + } + + //! 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; + } + +private: + //! Pointer to the single instance of the factory + static ReactionFactory* s_factory; + + //! default constructor, which is defined as private + ReactionFactory(); + + //! Mutex for use when calling the factory + static std::mutex reaction_mutex; + + //! 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. + */ +unique_ptr newReaction(const std::string& type); + +//! Create a new Reaction object for the reaction defined in `rxn_node` +/*! + * @param rxn_node XML node describing reaction. + */ +unique_ptr newReaction(const XML_Node& rxn_node); + +//! Create a new Reaction object using the specified parameters +/*! + * @param rxn_node AnyMap node describing reaction. + * @param kin kinetics manager + */ +unique_ptr newReaction(const AnyMap& rxn_node, + const Kinetics& kin); +} +#endif diff --git a/include/cantera/kinetics/ReactionRate.h b/include/cantera/kinetics/ReactionRate.h new file mode 100644 index 00000000000..52635465b84 --- /dev/null +++ b/include/cantera/kinetics/ReactionRate.h @@ -0,0 +1,203 @@ +/** + * @file ReactionRate.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 +// at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_REACTIONRATE_H +#define CT_REACTIONRATE_H + +#include "cantera/kinetics/RxnRates.h" +#include "cantera/kinetics/ReactionData.h" +#include "cantera/base/ctexceptions.h" + +namespace Cantera +{ + +class Func1; +class AnyMap; + + +//! Abstract base class for reaction rate definitions +/** + * Because this class has no template parameters, derived objects can be + * accessed via `shared_ptr`. For performance reasons + * it is essential that derived classes use the keyword `final` to + * de-virtualize `virtual` methods. + * + * 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: + virtual ~ReactionRateBase() {} + + //! Identifier of reaction type + virtual std::string type() const = 0; + + //! Evaluate reaction rate based on temperature + //! @param T temperature [K] + virtual double eval(double T) const = 0; + + //! 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; + + //! Evaluate reaction rate based on bulk phase + //! @param bulk object representing bulk phase + virtual double eval(const ThermoPhase& bulk) const = 0; + + //! Evaluate reaction rate derivative based on temperature + //! @param T temperature [K] + virtual double ddT(double T) const = 0; + + //! Evaluate reaction rate derivative based on temperature and pressure + //! @param T temperature [K] + //! @param P pressure [Pa] + virtual double ddT(double T, double P) const = 0; + + //! Evaluate reaction rate derivative based on bulk phase + //! @param bulk object representing bulk phase + virtual double ddT(const ThermoPhase& bulk) const = 0; +}; + + +//! Class template for reaction rate definitions with specialized DataType +/** + * 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 +{ +public: + ReactionRate() = default; + + //! 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) {} + + //! Evaluate reaction rate + //! @param shared_data data shared by all reactions of a given type + virtual double eval(const DataType& shared_data) const = 0; + + virtual double eval(double T) const override { + return eval(DataType(T)); + } + + virtual double eval(double T, double P) const override { + return eval(DataType(T, P)); + } + + virtual double eval(const ThermoPhase& bulk) const override { + return eval(DataType(bulk)); + } + + //! 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 { + throw CanteraError("ReactionRate::ddT", + "Not implemented by derived ReactionRate object."); + } + + virtual double ddT(double T) const override { + return ddT(DataType(T)); + } + + virtual double ddT(double T, double P) const override { + return ddT(DataType(T, P)); + } + + virtual double ddT(const ThermoPhase& bulk) const override { + return ddT(DataType(bulk)); + } +}; + + +//! Arrhenius reaction rate type depends only on temperature +/** + * Wrapped Arrhenius rate. + * + * @warning This class is an experimental part of the %Cantera API and + * may be changed or removed without notice. + */ +class ArrheniusRate final : public ReactionRate, public Arrhenius +{ +public: + ArrheniusRate(); + + 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 + ArrheniusRate(const AnyMap& node, const Units& rate_units); + + virtual std::string type() const override { return "ArrheniusRate"; } + + //! Update information specific to reaction + static bool uses_update() { return false; } + + virtual double eval(const ArrheniusData& shared_data) const override { + return updateRC(shared_data.m_logT, shared_data.m_recipT); + } + + virtual double ddT(const ArrheniusData& shared_data) const override { + return updateRC(shared_data.m_logT, shared_data.m_recipT) * + (m_b + m_E * shared_data.m_recipT) * shared_data.m_recipT; + } +}; + + +//! Custom reaction rate depending only on temperature +/** + * The rate expression is provided by a `Func1` object taking a single + * argument (temperature) and does not use a formalized parameterization. + * + * @warning This class is an experimental part of the %Cantera API and + * may be changed or removed without notice. + */ +class CustomFunc1Rate final : public ReactionRate +{ +public: + CustomFunc1Rate(); + + //! Constructor does nothing, as there is no formalized parameterization + //! @param node AnyMap object containing reaction rate specification + //! @param rate_units Description of units used for rate parameters + CustomFunc1Rate(const AnyMap& rate, const Units& rate_units) {} + + virtual std::string type() const override { return "custom-function"; } + + //! Update information specific to reaction + static bool uses_update() { return false; } + + //! Set custom rate + /** + * The call to the Func1 object takes a single argument (temperature) and + * does not depend on parameters handled in C++. + */ + void setRateFunction(shared_ptr f); + + virtual double eval(const CustomFunc1Data& shared_data) const override; + +protected: + shared_ptr m_ratefunc; +}; + +} + +#endif diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index c5e07de478d..7205eb62268 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -15,6 +15,9 @@ namespace Cantera { class Array2D; +class AnyValue; +class UnitSystem; +class Units; //! Arrhenius reaction rate type depends only on temperature /** @@ -38,6 +41,10 @@ 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 + void setParameters(const AnyValue& rate, + const UnitSystem& units, const Units& rate_units); + //! Update concentration-dependent parts of the rate coefficient. /*! * For this class, there are no concentration-dependent parts, so this diff --git a/include/cantera/kinetics/reaction_defs.h b/include/cantera/kinetics/reaction_defs.h index 8eaacdfbcea..7a1aedee97a 100644 --- a/include/cantera/kinetics/reaction_defs.h +++ b/include/cantera/kinetics/reaction_defs.h @@ -1,6 +1,8 @@ /** * @file reaction_defs.h * This file defines some constants used to specify reaction types. + * + * @deprecated To be removed after Cantera 2.6. */ // This file is part of Cantera. See License.txt in the top-level directory or @@ -57,6 +59,11 @@ 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 2a0f99b93f3..4b267fa4bc3 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -295,9 +295,12 @@ cdef extern from "cantera/thermo/SurfPhase.h": void getCoverages(double*) except +translate_exception -cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": +cdef extern from "cantera/kinetics/ReactionFactory.h" namespace "Cantera": + cdef shared_ptr[CxxReaction] CxxNewReaction "Cantera::newReaction" (string) except +translate_exception cdef shared_ptr[CxxReaction] CxxNewReaction "newReaction" (XML_Node&) except +translate_exception cdef shared_ptr[CxxReaction] CxxNewReaction "newReaction" (CxxAnyMap&, CxxKinetics&) except +translate_exception + +cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef vector[shared_ptr[CxxReaction]] CxxGetReactions "getReactions" (XML_Node&) except +translate_exception cdef vector[shared_ptr[CxxReaction]] CxxGetReactions "getReactions" (CxxAnyValue&, CxxKinetics&) except +translate_exception @@ -309,17 +312,33 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": double temperatureExponent() double activationEnergy_R() + cdef cppclass CxxReactionRateBase "Cantera::ReactionRateBase": + CxxReactionRateBase() + string type() + 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 + + cdef cppclass CxxArrheniusRate "Cantera::ArrheniusRate" (CxxReactionRateBase): + CxxArrheniusRate() + CxxArrheniusRate(double, double, double) + double preExponentialFactor() + double temperatureExponent() + double activationEnergy_R() + cdef cppclass CxxReaction "Cantera::Reaction": - # Note, this default constructor doesn't actually exist. The declaration - # is required by a Cython bug which should be resolved in Cython 0.22. CxxReaction() - CxxReaction(int) string reactantString() string productString() string equation() + string type() void validate() except +translate_exception - int reaction_type Composition reactants Composition products Composition orders @@ -329,6 +348,11 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cbool allow_nonreactant_orders cbool allow_negative_orders + cdef cppclass CxxReaction2 "Cantera::Reaction2" (CxxReaction): + CxxReaction2() + shared_ptr[CxxReactionRateBase] rate() + void setRate(shared_ptr[CxxReactionRateBase]) + cdef cppclass CxxElementaryReaction "Cantera::ElementaryReaction" (CxxReaction): CxxElementaryReaction() CxxArrhenius rate @@ -391,6 +415,13 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": 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 CxxCoverageDependency "Cantera::CoverageDependency": CxxCoverageDependency(double, double, double) double a @@ -431,6 +462,7 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera": shared_ptr[CxxReaction] reaction(size_t) except +translate_exception cbool isReversible(int) except +translate_exception int reactionType(int) except +translate_exception + string reactionTypeStr(int) except +translate_exception string reactionString(int) except +translate_exception string reactantString(int) except +translate_exception string productString(int) except +translate_exception @@ -1003,10 +1035,27 @@ cdef class ThermoPhase(_SolutionBase): cdef class InterfacePhase(ThermoPhase): cdef CxxSurfPhase* surf +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 Reaction: cdef shared_ptr[CxxReaction] _reaction cdef CxxReaction* reaction - cdef _assign(self, shared_ptr[CxxReaction] other) + @staticmethod + cdef wrap(shared_ptr[CxxReaction]) + +cdef class CustomReaction(Reaction): + cdef CustomRate _rate cdef class Arrhenius: cdef CxxArrhenius* rate @@ -1035,6 +1084,7 @@ cdef class Mixture: cpdef int element_index(self, element) except * cdef class Func1: + cdef shared_ptr[CxxFunc1] _func cdef CxxFunc1* func cdef object callable cdef object exception @@ -1184,7 +1234,6 @@ cdef np.ndarray get_transport_1d(Transport tran, transportMethod1d method) cdef np.ndarray get_transport_2d(Transport tran, transportMethod2d method) cdef CxxIdealGasPhase* getIdealGasPhase(ThermoPhase phase) except * cdef wrapSpeciesThermo(shared_ptr[CxxSpeciesThermo] spthermo) -cdef Reaction wrapReaction(shared_ptr[CxxReaction] reaction) cdef extern from "cantera/thermo/Elements.h" namespace "Cantera": double getElementWeight(string ename) except +translate_exception diff --git a/interfaces/cython/cantera/base.pyx b/interfaces/cython/cantera/base.pyx index d147980dcec..c92b008a0c4 100644 --- a/interfaces/cython/cantera/base.pyx +++ b/interfaces/cython/cantera/base.pyx @@ -223,7 +223,6 @@ cdef class _SolutionBase: for reaction in reactions: self.kinetics.addReaction(reaction._reaction) - def __getitem__(self, selection): copy = self.__class__(origin=self) if isinstance(selection, slice): diff --git a/interfaces/cython/cantera/examples/kinetics/custom_reactions.py b/interfaces/cython/cantera/examples/kinetics/custom_reactions.py new file mode 100644 index 00000000000..267d746569d --- /dev/null +++ b/interfaces/cython/cantera/examples/kinetics/custom_reactions.py @@ -0,0 +1,95 @@ +""" +An example demonstrating how to use custom reaction objects. + +For benchmark purposes, an ignition test is run to compare simulation times. + +Requires: cantera >= 2.6.0 +""" + +from timeit import default_timer +import numpy as np +from math import exp + +import cantera as ct + +mech = 'gri30.yaml' +fuel = 'CH4' +gas0 = ct.Solution(mech) + +species = gas0.species() +reactions = gas0.reactions() + +# construct custom reactions: replace 2nd reaction with equivalent custom reaction +custom_reactions = [r for r in reactions] +custom_reactions[2] = ct.CustomReaction( + equation='H2 + O <=> H + OH', + rate=lambda T: 38.7 * T**2.7 * exp(-3150.15428/T), + kinetics=gas0) + +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: + + 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) + +# construct test case - simulate ignition + +def ignition(gas): + # set up reactor + gas.TP = 1000., 5 * ct.one_atm + gas.set_equivalence_ratio(0.8, fuel, 'O2:1.0, N2:3.773') + r = ct.IdealGasReactor(gas) + net = ct.ReactorNet([r]) + net.rtol_sensitivity = 2.e-5 + + # time reactor integration + t1 = default_timer() + net.advance(.5) + t2 = default_timer() + + return 1000 * (t2 - t1) + +# output results + +repeat = 100 +print("Average time of {} simulation runs for '{}' " + "({})".format(repeat, mech, fuel)) + +sim0 = 0 +for i in range(repeat): + sim0 += ignition(gas0) +sim0 /= repeat +print('- Original mechanism: ' + '{0:.2f} ms (T_final={1:.2f})'.format(sim0, gas0.T)) + +sim1 = 0 +for i in range(repeat): + sim1 += ignition(gas1) +sim1 /= repeat +print('- One Python reaction: ' + '{0:.2f} ms (T_final={1:.2f}) ... ' + '{2:+.2f}%'.format(sim1, gas1.T, 100 * sim1 / sim0 - 100)) + +sim2 = 0 +for i in range(repeat): + sim2 += ignition(gas2) +sim2 /= repeat +print('- Alternative reactions: ' + '{0:.2f} ms (T_final={1:.2f}) ... ' + '{2:+.2f}%'.format(sim2, gas2.T, 100 * sim2 / sim0 - 100)) diff --git a/interfaces/cython/cantera/func1.pyx b/interfaces/cython/cantera/func1.pyx index 802537e3abc..9315576998d 100644 --- a/interfaces/cython/cantera/func1.pyx +++ b/interfaces/cython/cantera/func1.pyx @@ -86,10 +86,8 @@ cdef class Func1: cpdef void _set_callback(self, c) except *: self.callable = c - self.func = new CxxFunc1(func_callback, self) - - def __dealloc__(self): - del self.func + self._func.reset(new CxxFunc1(func_callback, self)) + self.func = self._func.get() def __call__(self, t): return self.func.eval(t) diff --git a/interfaces/cython/cantera/kinetics.pyx b/interfaces/cython/cantera/kinetics.pyx index d8a3f53e07a..866b9e5ce5e 100644 --- a/interfaces/cython/cantera/kinetics.pyx +++ b/interfaces/cython/cantera/kinetics.pyx @@ -105,7 +105,7 @@ cdef class Kinetics(_SolutionBase): ``i_reaction``. Changes to this object do not affect the `Kinetics` or `Solution` object until the `modify_reaction` function is called. """ - return wrapReaction(self.kinetics.reaction(i_reaction)) + return Reaction.wrap(self.kinetics.reaction(i_reaction)) def reactions(self): """ @@ -158,10 +158,23 @@ cdef class Kinetics(_SolutionBase): self.kinetics.setMultiplier(i_reaction, value) def reaction_type(self, int i_reaction): - """Type of reaction *i_reaction*.""" + """ + Type code of reaction *i_reaction*. + + .. deprecated:: 2.6 + """ + warnings.warn("Behavior changes after Cantera 2.6, " + "when function will return a type string. " + "Use method 'reaction_type_str' during " + "transition instead.", DeprecationWarning) self._check_reaction_index(i_reaction) return self.kinetics.reactionType(i_reaction) + def reaction_type_str(self, int i_reaction): + """Type of reaction *i_reaction*.""" + self._check_reaction_index(i_reaction) + return pystr(self.kinetics.reactionTypeStr(i_reaction)) + def reaction_equation(self, int i_reaction): """The equation for the specified reaction. See also `reaction_equations`.""" self._check_reaction_index(i_reaction) diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index 0f7bfa6de5b..ee2370c784d 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -1,14 +1,9 @@ # This file is part of Cantera. See License.txt in the top-level directory or # at https://cantera.org/license.txt for license and copyright information. -cdef extern from "cantera/kinetics/reaction_defs.h" namespace "Cantera": - cdef int ELEMENTARY_RXN - cdef int THREE_BODY_RXN - cdef int FALLOFF_RXN - cdef int PLOG_RXN - cdef int CHEBYSHEV_RXN - cdef int CHEMACT_RXN - cdef int INTERFACE_RXN + +# dictionary to store reaction classes +cdef dict _reaction_class_registry = {} cdef class Reaction: @@ -61,20 +56,43 @@ cdef class Reaction: R = ct.Reaction.fromCti('''reaction('O + H2 <=> H + OH', [(3.87e4, 'cm3/mol/s'), 2.7, (6260, 'cal/mol')])''') """ - reaction_type = 0 + reaction_type = "" def __cinit__(self, reactants='', products='', init=True, **kwargs): + if init: - self._reaction.reset(newReaction(self.reaction_type)) + self._reaction = CxxNewReaction(stringify((self.reaction_type))) self.reaction = self._reaction.get() if reactants: self.reactants = reactants if products: self.products = products - cdef _assign(self, shared_ptr[CxxReaction] other): - self._reaction = other - self.reaction = self._reaction.get() + @staticmethod + cdef wrap(shared_ptr[CxxReaction] reaction): + """ + Wrap a C++ Reaction object with a Python object of the correct derived type. + """ + # ensure all reaction types are registered + if not(_reaction_class_registry): + def register_subclasses(cls): + for c in cls.__subclasses__(): + _reaction_class_registry[getattr(c, 'reaction_type')] = c + register_subclasses(c) + + # update global reaction class registry + register_subclasses(Reaction) + + # identify class + rxn_type = pystr(reaction.get().type()) + cls = _reaction_class_registry.get(rxn_type, Reaction) + + # wrap C++ reaction + cdef Reaction R + R = cls(init=False) + R._reaction = reaction + R.reaction = R._reaction.get() + return R @staticmethod def fromCti(text): @@ -87,7 +105,7 @@ cdef class Reaction: """ cxx_reactions = CxxGetReactions(deref(CxxGetXmlFromString(stringify(text)))) assert cxx_reactions.size() == 1, cxx_reactions.size() - return wrapReaction(cxx_reactions[0]) + return Reaction.wrap(cxx_reactions[0]) @staticmethod def fromXml(text): @@ -99,7 +117,7 @@ cdef class Reaction: The XML input format is deprecated and will be removed in Cantera 3.0. """ cxx_reaction = CxxNewReaction(deref(CxxGetXmlFromString(stringify(text)))) - return wrapReaction(cxx_reaction) + return Reaction.wrap(cxx_reaction) @staticmethod def fromYaml(text, Kinetics kinetics): @@ -114,7 +132,7 @@ cdef class Reaction: """ cxx_reaction = CxxNewReaction(AnyMapFromYamlString(stringify(text)), deref(kinetics.kinetics)) - return wrapReaction(cxx_reaction) + return Reaction.wrap(cxx_reaction) @staticmethod def listFromFile(filename, Kinetics kinetics=None, section='reactions'): @@ -145,7 +163,7 @@ cdef class Reaction: deref(kinetics.kinetics)) else: cxx_reactions = CxxGetReactions(deref(CxxGetXmlFile(stringify(filename)))) - return [wrapReaction(r) for r in cxx_reactions] + return [Reaction.wrap(r) for r in cxx_reactions] @staticmethod def listFromXml(text): @@ -160,7 +178,7 @@ cdef class Reaction: The XML input format is deprecated and will be removed in Cantera 3.0. """ cxx_reactions = CxxGetReactions(deref(CxxGetXmlFromString(stringify(text)))) - return [wrapReaction(r) for r in cxx_reactions] + return [Reaction.wrap(r) for r in cxx_reactions] @staticmethod def listFromCti(text): @@ -175,7 +193,7 @@ cdef class Reaction: # Currently identical to listFromXml since get_XML_from_string is able # to distinguish between CTI and XML. cxx_reactions = CxxGetReactions(deref(CxxGetXmlFromString(stringify(text)))) - return [wrapReaction(r) for r in cxx_reactions] + return [Reaction.wrap(r) for r in cxx_reactions] @staticmethod def listFromYaml(text, Kinetics kinetics): @@ -186,7 +204,7 @@ cdef class Reaction: root = AnyMapFromYamlString(stringify(text)) cxx_reactions = CxxGetReactions(root[stringify("items")], deref(kinetics.kinetics)) - return [wrapReaction(r) for r in cxx_reactions] + return [Reaction.wrap(r) for r in cxx_reactions] property reactant_string: """ @@ -375,12 +393,155 @@ cdef copyArrhenius(CxxArrhenius* rate): 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 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:: + + rxn = ElementaryReaction(equation='H2 + O <=> H + OH', + rate={'A': 38.7, 'b': 2.7, 'Ea': 2.619184e+07}, + kinetics=gas) """ - reaction_type = ELEMENTARY_RXN + reaction_type = "elementary" + + 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, Arrhenius) 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, Arrhenius): + self.rate = rate property rate: """ Get/Set the `Arrhenius` rate coefficient for this reaction. """ @@ -409,7 +570,7 @@ cdef class ThreeBodyReaction(ElementaryReaction): A reaction with a non-reacting third body "M" that acts to add or remove energy from the reacting species. """ - reaction_type = THREE_BODY_RXN + reaction_type = "three-body" cdef CxxThreeBodyReaction* tbr(self): return self.reaction @@ -538,7 +699,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_RXN + reaction_type = "falloff" cdef CxxFalloffReaction* frxn(self): return self.reaction @@ -615,7 +776,7 @@ cdef class ChemicallyActivatedReaction(FalloffReaction): that the forward rate constant is written as being proportional to the low- pressure rate constant. """ - reaction_type = CHEMACT_RXN + reaction_type = "chemically-activated" cdef class PlogReaction(Reaction): @@ -623,7 +784,7 @@ cdef class PlogReaction(Reaction): A pressure-dependent reaction parameterized by logarithmically interpolating between Arrhenius rate expressions at various pressures. """ - reaction_type = PLOG_RXN + reaction_type = "pressure-dependent-Arrhenius" property rates: """ @@ -666,7 +827,7 @@ cdef class ChebyshevReaction(Reaction): A pressure-dependent reaction parameterized by a bivariate Chebyshev polynomial in temperature and pressure. """ - reaction_type = CHEBYSHEV_RXN + reaction_type = "Chebyshev" property Tmin: """ Minimum temperature [K] for the Chebyshev fit """ @@ -741,9 +902,106 @@ cdef class ChebyshevReaction(Reaction): return r.rate.updateRC(logT, recipT) +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_RXN + reaction_type = "interface" property coverage_deps: """ @@ -809,51 +1067,3 @@ cdef class InterfaceReaction(ElementaryReaction): def __set__(self, species): cdef CxxInterfaceReaction* r = self.reaction r.sticking_species = stringify(species) - - -cdef Reaction wrapReaction(shared_ptr[CxxReaction] reaction): - """ - Wrap a C++ Reaction object with a Python object of the correct derived type. - """ - cdef int reaction_type = reaction.get().reaction_type - - if reaction_type == ELEMENTARY_RXN: - R = ElementaryReaction(init=False) - elif reaction_type == THREE_BODY_RXN: - R = ThreeBodyReaction(init=False) - elif reaction_type == FALLOFF_RXN: - R = FalloffReaction(init=False) - elif reaction_type == CHEMACT_RXN: - R = ChemicallyActivatedReaction(init=False) - elif reaction_type == PLOG_RXN: - R = PlogReaction(init=False) - elif reaction_type == CHEBYSHEV_RXN: - R = ChebyshevReaction(init=False) - elif reaction_type == INTERFACE_RXN: - R = InterfaceReaction(init=False) - else: - R = Reaction(init=False) - - R._assign(reaction) - return R - -cdef CxxReaction* newReaction(int reaction_type): - """ - Create a new C++ Reaction object of the specified type - """ - if reaction_type == ELEMENTARY_RXN: - return new CxxElementaryReaction() - elif reaction_type == THREE_BODY_RXN: - return new CxxThreeBodyReaction() - elif reaction_type == FALLOFF_RXN: - return new CxxFalloffReaction() - elif reaction_type == CHEMACT_RXN: - return new CxxChemicallyActivatedReaction() - elif reaction_type == PLOG_RXN: - return new CxxPlogReaction() - elif reaction_type == CHEBYSHEV_RXN: - return new CxxChebyshevReaction() - elif reaction_type == INTERFACE_RXN: - return new CxxInterfaceReaction() - else: - return new CxxReaction(0) diff --git a/interfaces/cython/cantera/test/__init__.py b/interfaces/cython/cantera/test/__init__.py index 3eba767bd52..7e8905233c2 100644 --- a/interfaces/cython/cantera/test/__init__.py +++ b/interfaces/cython/cantera/test/__init__.py @@ -1,17 +1,18 @@ import os import cantera -from .test_thermo import * -from .test_purefluid import * +from .test_composite import * +from .test_convert import * from .test_equilibrium import * +from .test_func1 import * from .test_kinetics import * -from .test_transport import * from .test_mixture import * -from .test_func1 import * -from .test_reactor import * from .test_onedim import * -from .test_convert import * -from .test_composite import * +from .test_purefluid import * +from .test_reaction import * +from .test_reactor import * +from .test_thermo import * +from .test_transport import * cantera.add_directory(os.path.join(os.path.dirname(__file__), 'data')) cantera.add_directory(os.path.join(os.path.dirname(__file__), '..', 'examples', 'surface_chemistry')) diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 8b08fdff97b..7a21bce410b 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -50,14 +50,14 @@ def test_multiplier(self): self.assertArrayNear(0.5 * rev_rates0, rev_rates2) def test_reaction_type(self): - self.assertNear(self.phase.reaction_type(0), 2) # 3rd body - self.assertNear(self.phase.reaction_type(2), 1) # elementary - self.assertNear(self.phase.reaction_type(20), 4) # falloff + self.assertEqual(self.phase.reaction_type_str(0), "three-body") + self.assertEqual(self.phase.reaction_type_str(2), "elementary") + self.assertEqual(self.phase.reaction_type_str(20), "falloff") with self.assertRaisesRegex(ValueError, 'out of range'): - self.phase.reaction_type(33) + self.phase.reaction_type_str(33) with self.assertRaisesRegex(ValueError, 'out of range'): - self.phase.reaction_type(-2) + self.phase.reaction_type_str(-2) def test_reaction_equations(self): self.assertEqual(self.phase.n_reactions, diff --git a/interfaces/cython/cantera/test/test_reaction.py b/interfaces/cython/cantera/test/test_reaction.py new file mode 100644 index 00000000000..741c5ee7e44 --- /dev/null +++ b/interfaces/cython/cantera/test/test_reaction.py @@ -0,0 +1,140 @@ +from math import exp + +import cantera as ct +from . import utilities + + +class TestElementary(utilities.CanteraTest): + + _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" + + @classmethod + def setUpClass(cls): + utilities.CanteraTest.setUpClass() + cls.gas = ct.Solution('h2o2.xml') + 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.xml') + + def check_rxn(self, rxn): + 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) + + gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=self.species, reactions=[rxn]) + gas2.TPX = self.gas.TPX + self.check_sol(gas2) + + def check_sol(self, gas2): + ix = self._index + 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]) + + def test_from_parts(self): + orig = self.gas.reaction(self._index) + rxn = self._cls(orig.reactants, orig.products) + 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) + self.check_rxn(rxn) + + def test_from_rate(self): + rxn = self._cls(equation=self._equation, rate=self._rate_obj, kinetics=self.gas) + self.check_rxn(rxn) + + def test_add_rxn(self): + 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) + gas2.add_reaction(rxn) + self.check_sol(gas2) + + def test_wrong_rate(self): + with self.assertRaises(TypeError): + rxn = self._cls(equation=self._equation, rate=[], kinetics=self.gas) + + 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', + species=self.species, reactions=[rxn]) + gas2.TPX = self.gas.TPX + + self.assertNear(gas2.forward_rate_constants[0], 0.) + self.assertNear(gas2.net_rates_of_progress[0], 0.) + + def test_replace_rate(self): + rxn = self._cls(equation=self._equation, kinetics=self.gas) + rxn.rate = self._rate_obj + self.check_rxn(rxn) + + +class TestCustom(TestElementary): + + # 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" + + def setUp(self): + # 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) + + 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): + 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): + 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): + 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 7791ebe772c..f85dcba8180 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -120,6 +120,28 @@ bool BulkKinetics::addReaction(shared_ptr r) } else { m_irrev.push_back(nReactions()-1); } + + if (std::dynamic_pointer_cast(r) != nullptr) { + shared_ptr rate; + rate = std::dynamic_pointer_cast(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(); + + if (rate->type() == "ArrheniusRate") { + 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)); + } + } + + // Add reaction rate to evaluator + size_t index = m_bulk_types[rate->type()]; + m_bulk_rates[index]->add(nReactions() - 1, *rate); + } + return true; } @@ -128,6 +150,26 @@ void BulkKinetics::addElementaryReaction(ElementaryReaction& r) m_rates.install(nReactions()-1, r.rate); } +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(); + // Ensure that MultiBulkRates evaluator is available + if (m_bulk_types.find(rate->type()) != m_bulk_types.end()) { + throw CanteraError("BulkKinetics::modifyReaction", + "Evaluator not available for type '{}'.", rate->type()); + } + + // Replace reaction rate to evaluator + size_t index = m_bulk_types[rate->type()]; + m_bulk_rates[index]->replace(i, *rate); + } +} + void BulkKinetics::modifyElementaryReaction(size_t i, ElementaryReaction& rNew) { m_rates.replace(i, rNew.rate); diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 152e1f02538..ef008ab112f 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -23,10 +23,10 @@ GasKinetics::GasKinetics(ThermoPhase* thermo) : void GasKinetics::update_rates_T() { - doublereal T = thermo().temperature(); - doublereal P = thermo().pressure(); + double T = thermo().temperature(); + double P = thermo().pressure(); m_logStandConc = log(thermo().standardConcentration()); - doublereal logT = log(T); + double logT = log(T); if (T != m_temp) { if (!m_rfn.empty()) { @@ -45,6 +45,13 @@ void GasKinetics::update_rates_T() } if (T != m_temp || P != m_pres) { + + // loop over MultiBulkRates evaluators + for (auto& rates : m_bulk_rates) { + rates->update(thermo()); + rates->getRateConstants(thermo(), m_rfn.data()); + } + if (m_plog_rates.nReactions()) { m_plog_rates.update(T, logT, m_rfn.data()); m_ROP_ok = false; @@ -142,7 +149,7 @@ void GasKinetics::processFalloffReactions() m_falloffn.pr_to_falloff(pr.data(), falloff_work.data()); for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) { - if (reactionType(m_fallindx[i]) == FALLOFF_RXN) { + if (reactionTypeStr(m_fallindx[i]) == "falloff") { pr[i] *= m_rfn_high[i]; } else { // CHEMACT_RXN pr[i] *= m_rfn_low[i]; @@ -229,28 +236,26 @@ bool GasKinetics::addReaction(shared_ptr r) bool added = BulkKinetics::addReaction(r); if (!added) { return false; + } else if (std::dynamic_pointer_cast(r) != nullptr) { + // Rate object already added in BulkKinetics::addReaction + return true; } - switch (r->reaction_type) { - case ELEMENTARY_RXN: + if (r->type() == "elementary") { addElementaryReaction(dynamic_cast(*r)); - break; - case THREE_BODY_RXN: + } else if (r->type() == "three-body") { addThreeBodyReaction(dynamic_cast(*r)); - break; - case FALLOFF_RXN: - case CHEMACT_RXN: + } else if (r->type() == "falloff") { addFalloffReaction(dynamic_cast(*r)); - break; - case PLOG_RXN: + } else if (r->type() == "chemically-activated") { + addFalloffReaction(dynamic_cast(*r)); + } else if (r->type() == "pressure-dependent-Arrhenius") { addPlogReaction(dynamic_cast(*r)); - break; - case CHEBYSHEV_RXN: + } else if (r->type() == "Chebyshev") { addChebyshevReaction(dynamic_cast(*r)); - break; - default: + } else { throw CanteraError("GasKinetics::addReaction", - "Unknown reaction type specified: {}", r->reaction_type); + "Unknown reaction type specified: '{}'", r->type()); } return true; } @@ -286,7 +291,7 @@ void GasKinetics::addFalloffReaction(FalloffReaction& r) concm_falloff_values.resize(m_falloff_concm.workSize()); // install the falloff function calculator for this reaction - m_falloffn.install(nfall, r.reaction_type, r.falloff); + m_falloffn.install(nfall, r.type(), r.falloff); falloff_work.resize(m_falloffn.workSize()); } @@ -321,29 +326,29 @@ void GasKinetics::addChebyshevReaction(ChebyshevReaction& r) void GasKinetics::modifyReaction(size_t i, shared_ptr rNew) { - // operations common to all reaction types + // operations common to all bulk reaction types BulkKinetics::modifyReaction(i, rNew); - switch (rNew->reaction_type) { - case ELEMENTARY_RXN: + if (std::dynamic_pointer_cast(rNew) != nullptr) { + // Rate object already modified in BulkKinetics::modifyReaction + return; + } + + if (rNew->type() == "elementary") { modifyElementaryReaction(i, dynamic_cast(*rNew)); - break; - case THREE_BODY_RXN: + } else if (rNew->type() == "three-body") { modifyThreeBodyReaction(i, dynamic_cast(*rNew)); - break; - case FALLOFF_RXN: - case CHEMACT_RXN: + } else if (rNew->type() == "falloff") { + modifyFalloffReaction(i, dynamic_cast(*rNew)); + } else if (rNew->type() == "chemically-activated") { modifyFalloffReaction(i, dynamic_cast(*rNew)); - break; - case PLOG_RXN: + } else if (rNew->type() == "pressure-dependent-Arrhenius") { modifyPlogReaction(i, dynamic_cast(*rNew)); - break; - case CHEBYSHEV_RXN: + } else if (rNew->type() == "Chebyshev") { modifyChebyshevReaction(i, dynamic_cast(*rNew)); - break; - default: + } else { throw CanteraError("GasKinetics::modifyReaction", - "Unknown reaction type specified: {}", rNew->reaction_type); + "Unknown reaction type specified: '{}'", rNew->type()); } // invalidate all cached data diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 42ab44e0386..9f883583a51 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -14,6 +14,7 @@ #include "cantera/thermo/ThermoPhase.h" #include "cantera/base/stringUtils.h" #include "cantera/base/utilities.h" +#include "cantera/base/global.h" #include using namespace std; @@ -115,7 +116,7 @@ std::pair Kinetics::checkDuplicates(bool throw_err) const unmatched_duplicates.erase(i); unmatched_duplicates.erase(m); continue; - } else if (R.reaction_type != other.reaction_type) { + } else if (R.type() != other.type()) { continue; // different reaction types } doublereal c = checkDuplicateStoich(net_stoich[i], net_stoich[m]); @@ -123,8 +124,8 @@ std::pair Kinetics::checkDuplicates(bool throw_err) const continue; // stoichiometries differ (not by a multiple) } else if (c < 0.0 && !R.reversible && !other.reversible) { continue; // irreversible reactions in opposite directions - } else if (R.reaction_type == FALLOFF_RXN || - R.reaction_type == CHEMACT_RXN) { + } else if (R.type() == "falloff" || + R.type() == "chemically-activated") { ThirdBody& tb1 = dynamic_cast(R).third_body; ThirdBody& tb2 = dynamic_cast(other).third_body; bool thirdBodyOk = true; @@ -140,7 +141,7 @@ std::pair Kinetics::checkDuplicates(bool throw_err) const if (thirdBodyOk) { continue; // No overlap in third body efficiencies } - } else if (R.reaction_type == THREE_BODY_RXN) { + } else if (R.type() == "three-body") { ThirdBody& tb1 = dynamic_cast(R).third_body; ThirdBody& tb2 = dynamic_cast(other).third_body; bool thirdBodyOk = true; @@ -370,9 +371,16 @@ double Kinetics::productStoichCoeff(size_t kSpec, size_t irxn) const } int Kinetics::reactionType(size_t i) const { + warn_deprecated("Kinetics::reactionType()", + "To be changed after Cantera 2.6. " + "Return string instead of magic number; use " + "Kinetics::reactionTypeStr during transition."); return m_reactions[i]->reaction_type; } +std::string Kinetics::reactionTypeStr(size_t i) const { + return m_reactions[i]->type(); +} std::string Kinetics::reactionString(size_t i) const { @@ -613,10 +621,10 @@ void Kinetics::modifyReaction(size_t i, shared_ptr rNew) { checkReactionIndex(i); shared_ptr& rOld = m_reactions[i]; - if (rNew->reaction_type != rOld->reaction_type) { + if (rNew->type() != rOld->type()) { throw CanteraError("Kinetics::modifyReaction", "Reaction types are different: {} != {}.", - rOld->reaction_type, rNew->reaction_type); + rOld->type(), rNew->type()); } if (rNew->reactants != rOld->reactants) { diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index c5f82e8c595..23b5047ba70 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -6,6 +6,7 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/ReactionFactory.h" #include "cantera/kinetics/FalloffFactory.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/thermo/ThermoPhase.h" @@ -24,13 +25,40 @@ namespace ba = boost::algorithm; namespace Cantera { +Reaction::Reaction() + : reaction_type(NONE) + , reversible(true) + , duplicate(false) + , allow_nonreactant_orders(false) + , allow_negative_orders(false) + , m_valid(true) +{ +} + +Reaction::Reaction(const Composition& reactants_, + const Composition& products_) + : reaction_type(NONE) + , reactants(reactants_) + , products(products_) + , reversible(true) + , duplicate(false) + , allow_nonreactant_orders(false) + , allow_negative_orders(false) + , m_valid(true) +{ +} + Reaction::Reaction(int type) : reaction_type(type) , reversible(true) , duplicate(false) , allow_nonreactant_orders(false) , allow_negative_orders(false) + , m_valid(true) { + warn_deprecated("Reaction::Reaction()", + "To be removed after Cantera 2.6. Use constructor without parameter " + "'type' instead."); } Reaction::Reaction(int type, const Composition& reactants_, @@ -42,7 +70,11 @@ Reaction::Reaction(int type, const Composition& reactants_, , duplicate(false) , allow_nonreactant_orders(false) , allow_negative_orders(false) + , m_valid(true) { + warn_deprecated("Reaction::Reaction()", + "To be removed after Cantera 2.6. Use constructor without parameter " + "'type' instead."); } void Reaction::validate() @@ -108,16 +140,18 @@ std::string Reaction::equation() const ElementaryReaction::ElementaryReaction(const Composition& reactants_, const Composition products_, const Arrhenius& rate_) - : Reaction(ELEMENTARY_RXN, reactants_, products_) + : Reaction(reactants_, products_) , rate(rate_) , allow_negative_pre_exponential_factor(false) { + reaction_type = ELEMENTARY_RXN; } ElementaryReaction::ElementaryReaction() - : Reaction(ELEMENTARY_RXN) + : Reaction() , allow_negative_pre_exponential_factor(false) { + reaction_type = ELEMENTARY_RXN; } void ElementaryReaction::validate() @@ -165,22 +199,24 @@ std::string ThreeBodyReaction::productString() const { } FalloffReaction::FalloffReaction() - : Reaction(FALLOFF_RXN) + : Reaction() , falloff(new Falloff()) , allow_negative_pre_exponential_factor(false) { + reaction_type = FALLOFF_RXN; } FalloffReaction::FalloffReaction( const Composition& reactants_, const Composition& products_, const Arrhenius& low_rate_, const Arrhenius& high_rate_, const ThirdBody& tbody) - : Reaction(FALLOFF_RXN, reactants_, products_) + : Reaction(reactants_, products_) , low_rate(low_rate_) , high_rate(high_rate_) , third_body(tbody) , falloff(new Falloff()) { + reaction_type = FALLOFF_RXN; } std::string FalloffReaction::reactantString() const { @@ -233,28 +269,87 @@ ChemicallyActivatedReaction::ChemicallyActivatedReaction( } PlogReaction::PlogReaction() - : Reaction(PLOG_RXN) + : Reaction() { + reaction_type = PLOG_RXN; } PlogReaction::PlogReaction(const Composition& reactants_, const Composition& products_, const Plog& rate_) - : Reaction(PLOG_RXN, reactants_, products_) + : Reaction(reactants_, products_) , rate(rate_) { + reaction_type = PLOG_RXN; } ChebyshevReaction::ChebyshevReaction() - : Reaction(CHEBYSHEV_RXN) + : Reaction() { + reaction_type = CHEBYSHEV_RXN; } ChebyshevReaction::ChebyshevReaction(const Composition& reactants_, const Composition& products_, const ChebyshevRate& rate_) - : Reaction(CHEBYSHEV_RXN, reactants_, products_) + : 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); + + input = node; +} + +CustomFunc1Reaction::CustomFunc1Reaction() + : Reaction2() +{ + reaction_type = CUSTOMPY_RXN; +} + +void CustomFunc1Reaction::setParameters(const AnyMap& node, const Kinetics& kin) +{ + Reaction2::setParameters(node, kin); + Units rate_units; // @todo Not needed once `rate_units` is implemented. + 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); + + // @todo Rate units will become available as `rate_units` after + // serialization is fully implemented. + Units rate_units = rateCoeffUnits(*this, kin); + setRate( + std::shared_ptr(new ArrheniusRate(node, rate_units))); + allow_negative_pre_exponential_factor = node.getBool("negative-A", false); } InterfaceReaction::InterfaceReaction() @@ -298,13 +393,13 @@ Arrhenius readArrhenius(const XML_Node& arrhenius_node) } Units rateCoeffUnits(const Reaction& R, const Kinetics& kin, - int pressure_dependence=0) + int pressure_dependence) { - if (R.reaction_type == INVALID_RXN) { + if (!R.valid()) { // If a reaction is invalid because of missing species in the Kinetics // object, determining the units of the rate coefficient is impossible. return Units(); - } else if (R.reaction_type == INTERFACE_RXN + } else if (R.type() == "interface" && dynamic_cast(R).is_sticking_coefficient) { // Sticking coefficients are dimensionless return Units(); @@ -339,20 +434,10 @@ Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate, const Kinetics& kin, const UnitSystem& units, int pressure_dependence=0) { - double A, b, Ta; - Units rc_units = rateCoeffUnits(R, kin, pressure_dependence); - if (rate.is()) { - auto& rate_map = rate.as(); - A = units.convert(rate_map["A"], rc_units); - b = rate_map["b"].asDouble(); - Ta = units.convertActivationEnergy(rate_map["Ea"], "K"); - } else { - auto& rate_vec = rate.asVector(3); - A = units.convert(rate_vec[0], rc_units); - b = rate_vec[1].asDouble(); - Ta = units.convertActivationEnergy(rate_vec[2], "K"); - } - return Arrhenius(A, b, Ta); + Units rate_units = rateCoeffUnits(R, kin, pressure_dependence); + Arrhenius arr; + arr.setParameters(rate, units, rate_units); + return arr; } //! Parse falloff parameters, given a rateCoeff node @@ -508,7 +593,7 @@ void parseReactionEquation(Reaction& R, const AnyValue& equation, } if (kin.kineticsSpeciesIndex(species) == npos && stoich != -1 && species != "M") { - R.reaction_type = INVALID_RXN; + R.setValid(false); } if (reactants) { @@ -540,7 +625,7 @@ void setupReaction(Reaction& R, const AnyMap& node, const Kinetics& kin) for (const auto& order : node["orders"].asMap()) { R.orders[order.first] = order.second; if (kin.kineticsSpeciesIndex(order.first) == npos) { - R.reaction_type = INVALID_RXN; + R.setValid(false); } } } @@ -906,141 +991,6 @@ void setupElectrochemicalReaction(ElectrochemicalReaction& R, "exchange-current-density-formulation", false); } -bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin) -{ - vector_fp e_counter(kin.nPhases(), 0.0); - - // Find the number of electrons in the products for each phase - for (const auto& sp : R.products) { - size_t kkin = kin.kineticsSpeciesIndex(sp.first); - size_t i = kin.speciesPhaseIndex(kkin); - size_t kphase = kin.thermo(i).speciesIndex(sp.first); - e_counter[i] += sp.second * kin.thermo(i).charge(kphase); - } - - // Subtract the number of electrons in the reactants for each phase - for (const auto& sp : R.reactants) { - size_t kkin = kin.kineticsSpeciesIndex(sp.first); - size_t i = kin.speciesPhaseIndex(kkin); - size_t kphase = kin.thermo(i).speciesIndex(sp.first); - e_counter[i] -= sp.second * kin.thermo(i).charge(kphase); - } - - // If the electrons change phases then the reaction is electrochemical - for (double delta_e : e_counter) { - if (std::abs(delta_e) > 1e-4) { - return true; - } - } - return false; -} - -shared_ptr newReaction(const XML_Node& rxn_node) -{ - std::string type = toLowerCopy(rxn_node["type"]); - - // Modify the reaction type for interface reactions which contain - // electrochemical reaction data - if (rxn_node.child("rateCoeff").hasChild("electrochem") - && (type == "edge" || type == "surface")) { - type = "electrochemical"; - } - - // Create a new Reaction object of the appropriate type - if (type == "elementary" || type == "arrhenius" || type == "") { - auto R = make_shared(); - setupElementaryReaction(*R, rxn_node); - return R; - } else if (type == "threebody" || type == "three_body") { - auto R = make_shared(); - setupThreeBodyReaction(*R, rxn_node); - return R; - } else if (type == "falloff") { - auto R = make_shared(); - setupFalloffReaction(*R, rxn_node); - return R; - } else if (type == "chemact" || type == "chemically_activated") { - auto R = make_shared(); - setupChemicallyActivatedReaction(*R, rxn_node); - return R; - } else if (type == "plog" || type == "pdep_arrhenius") { - auto R = make_shared(); - setupPlogReaction(*R, rxn_node); - return R; - } else if (type == "chebyshev") { - auto R = make_shared(); - setupChebyshevReaction(*R, rxn_node); - return R; - } else if (type == "interface" || type == "surface" || type == "edge" || - type == "global") { - auto R = make_shared(); - setupInterfaceReaction(*R, rxn_node); - return R; - } else if (type == "electrochemical" || - type == "butlervolmer_noactivitycoeffs" || - type == "butlervolmer" || - type == "surfaceaffinity") { - auto R = make_shared(); - setupElectrochemicalReaction(*R, rxn_node); - return R; - } else { - throw CanteraError("newReaction", - "Unknown reaction type '" + rxn_node["type"] + "'"); - } -} - -unique_ptr newReaction(const AnyMap& node, const Kinetics& kin) -{ - std::string type = "elementary"; - if (node.hasKey("type")) { - type = node["type"].asString(); - } - - if (kin.thermo(kin.reactionPhaseIndex()).nDim() < 3) { - // See if this is an electrochemical reaction - Reaction testReaction(0); - parseReactionEquation(testReaction, node["equation"], kin); - if (isElectrochemicalReaction(testReaction, kin)) { - unique_ptr R(new ElectrochemicalReaction()); - setupElectrochemicalReaction(*R, node, kin); - return unique_ptr(move(R)); - } else { - unique_ptr R(new InterfaceReaction()); - setupInterfaceReaction(*R, node, kin); - return unique_ptr(move(R)); - } - } - - if (type == "elementary") { - unique_ptr R(new ElementaryReaction()); - setupElementaryReaction(*R, node, kin); - return unique_ptr(move(R)); - } else if (type == "three-body") { - unique_ptr R(new ThreeBodyReaction()); - setupThreeBodyReaction(*R, node, kin); - return unique_ptr(move(R)); - } else if (type == "falloff") { - unique_ptr R(new FalloffReaction()); - setupFalloffReaction(*R, node, kin); - return unique_ptr(move(R)); - } else if (type == "chemically-activated") { - unique_ptr R(new ChemicallyActivatedReaction()); - setupFalloffReaction(*R, node, kin); - return unique_ptr(move(R)); - } else if (type == "pressure-dependent-Arrhenius") { - unique_ptr R(new PlogReaction()); - setupPlogReaction(*R, node, kin); - return unique_ptr(move(R)); - } else if (type == "Chebyshev") { - unique_ptr R(new ChebyshevReaction()); - setupChebyshevReaction(*R, node, kin); - return unique_ptr(move(R)); - } else { - throw InputFileError("newReaction", node["type"], - "Unknown reaction type '{}'", type); - } -} - std::vector > getReactions(const XML_Node& node) { std::vector > all_reactions; @@ -1056,7 +1006,7 @@ std::vector> getReactions(const AnyValue& items, std::vector> all_reactions; for (const auto& node : items.asVector()) { shared_ptr R(newReaction(node, kinetics)); - if (R->reaction_type != INVALID_RXN) { + if (R->valid()) { all_reactions.emplace_back(R); } else if (!kinetics.skipUndeclaredSpecies()) { throw InputFileError("getReactions", node, diff --git a/src/kinetics/ReactionData.cpp b/src/kinetics/ReactionData.cpp new file mode 100644 index 00000000000..07041e33a21 --- /dev/null +++ b/src/kinetics/ReactionData.cpp @@ -0,0 +1,20 @@ +//! @file ReactionData.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/kinetics/ReactionData.h" +#include "cantera/thermo/ThermoPhase.h" + +namespace Cantera +{ + +void ArrheniusData::update(const ThermoPhase& bulk) { + update(bulk.temperature()); +} + +void CustomFunc1Data::update(const ThermoPhase& bulk) { + m_temperature = bulk.temperature(); +} + +} diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp new file mode 100644 index 00000000000..fa741741f67 --- /dev/null +++ b/src/kinetics/ReactionFactory.cpp @@ -0,0 +1,241 @@ + /** + * @file ReactionFactory.cpp + */ + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/ReactionFactory.h" +#include "cantera/kinetics/Kinetics.h" +#include "cantera/base/ctml.h" +#include "cantera/base/AnyMap.h" +#include "cantera/base/stringUtils.h" + +namespace Cantera +{ + +ReactionFactory* ReactionFactory::s_factory = 0; +std::mutex ReactionFactory::reaction_mutex; + +ReactionFactory::ReactionFactory() +{ + // register elementary reactions + reg("elementary", []() { return new ElementaryReaction(); }); + 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 three-body reactions + reg("three-body", []() { return new ThreeBodyReaction(); }); + 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 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); + }); + + // register falloff reactions + reg("chemically-activated", []() { return new ChemicallyActivatedReaction(); }); + 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(); }); + 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 Chebyshev reactions + reg("Chebyshev", []() { return new ChebyshevReaction(); }); + 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); + }); + + // 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); + }); + + // register interface reactions + reg("interface", []() { return new InterfaceReaction(); }); + 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); + }); +} + +bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin) +{ + vector_fp e_counter(kin.nPhases(), 0.0); + + // Find the number of electrons in the products for each phase + for (const auto& sp : R.products) { + size_t kkin = kin.kineticsSpeciesIndex(sp.first); + size_t i = kin.speciesPhaseIndex(kkin); + size_t kphase = kin.thermo(i).speciesIndex(sp.first); + e_counter[i] += sp.second * kin.thermo(i).charge(kphase); + } + + // Subtract the number of electrons in the reactants for each phase + for (const auto& sp : R.reactants) { + size_t kkin = kin.kineticsSpeciesIndex(sp.first); + size_t i = kin.speciesPhaseIndex(kkin); + size_t kphase = kin.thermo(i).speciesIndex(sp.first); + e_counter[i] -= sp.second * kin.thermo(i).charge(kphase); + } + + // If the electrons change phases then the reaction is electrochemical + for (double delta_e : e_counter) { + if (std::abs(delta_e) > 1e-4) { + return true; + } + } + return false; +} + +unique_ptr newReaction(const std::string& type) +{ + unique_ptr R(ReactionFactory::factory()->create(type)); + return R; +} + +unique_ptr newReaction(const XML_Node& rxn_node) +{ + std::string type = toLowerCopy(rxn_node["type"]); + + // Modify the reaction type for interface reactions which contain + // electrochemical reaction data + if (rxn_node.child("rateCoeff").hasChild("electrochem") + && (type == "edge" || type == "surface")) { + type = "electrochemical"; + } + + Reaction* R; + try { + R = ReactionFactory::factory()->create(type); + } catch (CanteraError& err) { + throw CanteraError("newReaction", + "Unknown reaction type '" + rxn_node["type"] + "'"); + } + if (type != "electrochemical") { + type = R->type(); + } + ReactionFactory::factory()->setup_XML(type, R, rxn_node); + + return unique_ptr(R); +} + +unique_ptr newReaction(const AnyMap& rxn_node, + const Kinetics& kin) +{ + std::string type = "elementary"; + if (rxn_node.hasKey("type")) { + type = rxn_node["type"].asString(); + } + + if (kin.thermo().nDim() < 3) { + // See if this is an electrochemical reaction: type of + // receiving reaction object is unimportant in this case + ElementaryReaction testReaction; + parseReactionEquation(testReaction, rxn_node["equation"], kin); + if (isElectrochemicalReaction(testReaction, kin)) { + type = "electrochemical"; + } else { + type = "interface"; + } + } + + Reaction* R; + try { + R = ReactionFactory::factory()->create(type); + } catch (CanteraError& err) { + throw InputFileError("ReactionFactory::newReaction", rxn_node["type"], + "Unknown reaction type '{}'", type); + } + ReactionFactory::factory()->setup_AnyMap(type, R, rxn_node, kin); + + return unique_ptr(R); +} + +} diff --git a/src/kinetics/ReactionPath.cpp b/src/kinetics/ReactionPath.cpp index 7be482fdb96..9672e1cc3f1 100644 --- a/src/kinetics/ReactionPath.cpp +++ b/src/kinetics/ReactionPath.cpp @@ -680,9 +680,9 @@ string reactionLabel(size_t i, size_t kr, size_t nr, label += " + "+ s.kineticsSpeciesName(slist[j]); } } - if (s.reactionType(i) == THREE_BODY_RXN) { + if (s.reactionTypeStr(i) == "three-body") { label += " + M "; - } else if (s.reactionType(i) == FALLOFF_RXN) { + } else if (s.reactionTypeStr(i) == "falloff") { label += " (+ M)"; } return label; @@ -735,9 +735,9 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, revlabel += " + "+ s.kineticsSpeciesName(m_prod[i][j]); } } - if (s.reactionType(i) == THREE_BODY_RXN) { + if (s.reactionTypeStr(i) == "three-body") { revlabel += " + M "; - } else if (s.reactionType(i) == FALLOFF_RXN) { + } else if (s.reactionTypeStr(i) == "falloff") { revlabel += " (+ M)"; } diff --git a/src/kinetics/ReactionRate.cpp b/src/kinetics/ReactionRate.cpp new file mode 100644 index 00000000000..29485704c97 --- /dev/null +++ b/src/kinetics/ReactionRate.cpp @@ -0,0 +1,35 @@ +//! @file ReactionRate.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/kinetics/ReactionRate.h" +#include "cantera/numerics/Func1.h" +#include "cantera/base/AnyMap.h" + +namespace Cantera +{ + +ArrheniusRate::ArrheniusRate(double A, double b, double E) + : Arrhenius(A, b, E) { +} + +ArrheniusRate::ArrheniusRate(const AnyMap& node, const Units& rate_units) { + setParameters(node["rate-constant"], node.units(), rate_units); +} + +CustomFunc1Rate::CustomFunc1Rate() : m_ratefunc(0) {} + +void CustomFunc1Rate::setRateFunction(shared_ptr f) { + m_ratefunc = f; +} + +double CustomFunc1Rate::eval(const CustomFunc1Data& shared_data) const { + if (m_ratefunc) { + return m_ratefunc->eval(shared_data.m_temperature); + } + throw CanteraError("CustomFunc1Rate::eval", + "Custom rate function is not initialized."); +} + +} diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 269378b7729..0b722881587 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -5,6 +5,7 @@ #include "cantera/kinetics/RxnRates.h" #include "cantera/base/Array.h" +#include "cantera/base/AnyMap.h" namespace Cantera { @@ -28,6 +29,28 @@ Arrhenius::Arrhenius(doublereal A, doublereal b, doublereal E) } } +void Arrhenius::setParameters(const AnyValue& rate, + const UnitSystem& units, const Units& rate_units) +{ + if (rate.is()) { + auto& rate_map = rate.as(); + m_A = units.convert(rate_map["A"], rate_units); + m_b = rate_map["b"].asDouble(); + m_E = units.convertActivationEnergy(rate_map["Ea"], "K"); + } else { + auto& rate_vec = rate.asVector(3); + m_A = units.convert(rate_vec[0], rate_units); + m_b = rate_vec[1].asDouble(); + m_E = units.convertActivationEnergy(rate_vec[2], "K"); + } + + if (m_A <= 0.0) { + m_logA = -1.0E300; + } else { + m_logA = std::log(m_A); + } +} + SurfaceArrhenius::SurfaceArrhenius() : m_b(0.0) , m_E(0.0) diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index 296eb6f35dc..3164ead2b18 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -16,6 +16,7 @@ #include "cantera/kinetics/importKinetics.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/ReactionFactory.h" #include "cantera/base/stringUtils.h" #include "cantera/base/ctml.h" #include "cantera/base/yaml.h" diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 352fb5ff6d8..5ce804a7d71 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -4,6 +4,7 @@ #include "cantera/kinetics/GasKinetics.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/kinetics/KineticsFactory.h" +#include "cantera/kinetics/ReactionFactory.h" #include "cantera/thermo/ThermoFactory.h" using namespace Cantera; @@ -19,7 +20,7 @@ TEST(Reaction, ElementaryFromYaml) auto R = newReaction(rxn, *(sol->kinetics())); EXPECT_EQ(R->reactants.at("NO"), 1); EXPECT_EQ(R->products.at("N2"), 1); - EXPECT_EQ(R->reaction_type, ELEMENTARY_RXN); + EXPECT_EQ(R->type(), "elementary"); auto ER = dynamic_cast(*R); EXPECT_DOUBLE_EQ(ER.rate.preExponentialFactor(), -2.7e10); @@ -39,6 +40,7 @@ TEST(Reaction, ThreeBodyFromYaml1) auto R = newReaction(rxn, *(sol->kinetics())); EXPECT_EQ(R->reactants.count("M"), (size_t) 0); + EXPECT_EQ(R->type(), "three-body"); auto TBR = dynamic_cast(*R); EXPECT_DOUBLE_EQ(TBR.rate.preExponentialFactor(), 1.2e11);