diff --git a/include/cantera/base/ct_defs.h b/include/cantera/base/ct_defs.h index 8ba66218d69..7f7de6db776 100644 --- a/include/cantera/base/ct_defs.h +++ b/include/cantera/base/ct_defs.h @@ -106,6 +106,7 @@ const double ElectronMass = 9.1093837015e-31; //! @{ //! Reduced Planck constant \f$ \hbar \f$ [m2-kg/s] +//! @deprecated Unused. To be removed after Cantera 2.6. const double Planck_bar = Planck / (2 * Pi); //! Universal Gas Constant \f$ R_u \f$ [J/kmol/K] @@ -117,10 +118,11 @@ const double logGasConstant = std::log(GasConstant); const double GasConst_cal_mol_K = GasConstant / 4184.0; //! log(k_b/h) +//! @deprecated Unused. To be removed after Cantera 2.6. const double logBoltz_Planck = std::log(Boltzmann / Planck); //! Stefan-Boltzmann constant \f$ \sigma \f$ [W/m2/K4] -const double StefanBoltz = Pi * Pi * std::pow(Boltzmann, 4.0) / (60.0 * std::pow(Planck_bar, 3.0) * lightSpeed * lightSpeed); // 5.670374419e-8 +const double StefanBoltz = 2.0 * std::pow(Pi, 5) * std::pow(Boltzmann, 4) / (15.0 * std::pow(Planck, 3) * lightSpeed * lightSpeed); // 5.670374419e-8 //! Faraday constant \f$ F \f$ [C/kmol] const double Faraday = ElectronCharge * Avogadro; diff --git a/include/cantera/base/global.h b/include/cantera/base/global.h index 5e8da9b66a7..87a27a3386d 100644 --- a/include/cantera/base/global.h +++ b/include/cantera/base/global.h @@ -263,7 +263,7 @@ void suppress_thermo_warnings(bool suppress=true); //! @copydoc Application::thermo_warnings_suppressed bool thermo_warnings_suppressed(); -//! @copydoc Application::suppress_user_warnings +//! @copydoc Application::suppress_warnings void suppress_warnings(); //! @copydoc Application::warnings_suppressed diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index eb74573aa2f..fa59c986119 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -21,6 +21,18 @@ namespace Cantera class AnyValue; class AnyMap; +//! Data container holding shared data specific to ArrheniusRate +/** + * The data container `ArrheniusData` holds precalculated data common to + * all `ArrheniusRate` objects. + */ +struct ArrheniusData : public ReactionData +{ + virtual bool update(const ThermoPhase& phase, const Kinetics& kin); + using ReactionData::update; +}; + + /** * @defgroup arrheniusGroup Arrhenius-type Parameterizations * @@ -58,6 +70,12 @@ class ArrheniusBase : public ReactionRate setRateParameters(rate, units, rate_units); } + explicit ArrheniusBase(const AnyMap& node, const UnitStack& rate_units={}) + : ArrheniusBase() + { + setParameters(node, rate_units); + } + //! Perform object setup based on AnyValue node information /*! * @param rate AnyValue containing rate information @@ -71,9 +89,16 @@ class ArrheniusBase : public ReactionRate //! Return parameters void getRateParameters(AnyMap& node) const; + virtual void setParameters( + const AnyMap& node, const UnitStack& rate_units) override; + + virtual void getParameters(AnyMap& node) const override; + //! Check rate expression virtual void check(const std::string& equation, const AnyMap& node) override; + virtual void validate(const std::string& equation, const Kinetics& kin) override; + //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending //! on the reaction order) /*! @@ -149,17 +174,6 @@ class ArrheniusBase : public ReactionRate Units m_rate_units; //!< Reaction rate units }; -//! Data container holding shared data specific to ArrheniusRate -/** - * The data container `ArrheniusData` holds precalculated data common to - * all `ArrheniusRate` objects. - */ -struct ArrheniusData : public ReactionData -{ - virtual bool update(const ThermoPhase& phase, const Kinetics& kin); - using ReactionData::update; -}; - //! Arrhenius reaction rate type depends only on temperature /*! * A reaction rate coefficient of the following form. @@ -170,16 +184,19 @@ struct ArrheniusData : public ReactionData * * @ingroup arrheniusGroup * - * @todo supersedes Arrhenius2 and will replace Arrhenius(2) after Cantera 2.6, - * The class should be renamed to Arrhenius after removal of Arrhenius2. The new + * @todo supersedes Arrhenius2 and will replace Arrhenius after Cantera 2.6. The new * behavior can be forced in self-compiled Cantera installations by defining * CT_NO_LEGACY_REACTIONS_26 via the 'no_legacy_reactions' option in SCons. */ -class Arrhenius3 : public ArrheniusBase +class ArrheniusRate : public ArrheniusBase { public: using ArrheniusBase::ArrheniusBase; // inherit constructors + unique_ptr newMultiRate() const override { + return unique_ptr(new MultiRate); + } + virtual const std::string type() const override { return "Arrhenius"; } @@ -212,54 +229,6 @@ class Arrhenius3 : public ArrheniusBase } }; - -//! A class template for bulk phase reaction rate specifications -template -class BulkRate : public RateType -{ -public: - BulkRate() = default; - using RateType::RateType; // inherit constructors - - //! Constructor based on AnyMap content - BulkRate(const AnyMap& node, const UnitStack& rate_units={}) { - setParameters(node, rate_units); - } - - unique_ptr newMultiRate() const override { - return unique_ptr( - new MultiRate, DataType>); - } - - virtual void setParameters( - const AnyMap& node, const UnitStack& rate_units) override - { - RateType::m_negativeA_ok = node.getBool("negative-A", false); - if (!node.hasKey("rate-constant")) { - RateType::setRateParameters(AnyValue(), node.units(), rate_units); - return; - } - RateType::setRateParameters(node["rate-constant"], node.units(), rate_units); - } - - virtual void getParameters(AnyMap& node) const override { - if (RateType::m_negativeA_ok) { - node["negative-A"] = true; - } - AnyMap rateNode; - RateType::getRateParameters(rateNode); - if (!rateNode.empty()) { - // RateType object is configured - node["rate-constant"] = std::move(rateNode); - } - if (RateType::type() != "Arrhenius") { - node["type"] = RateType::type(); - } - } -}; - -typedef BulkRate ArrheniusRate; - } #endif diff --git a/include/cantera/kinetics/BlowersMaselRate.h b/include/cantera/kinetics/BlowersMaselRate.h index 23cc818e270..c4ea6f4ad80 100644 --- a/include/cantera/kinetics/BlowersMaselRate.h +++ b/include/cantera/kinetics/BlowersMaselRate.h @@ -67,11 +67,11 @@ struct BlowersMaselData : public ReactionData * * @ingroup arrheniusGroup */ -class BlowersMasel : public ArrheniusBase +class BlowersMaselRate : public ArrheniusBase { public: //! Default constructor. - BlowersMasel(); + BlowersMaselRate(); //! Constructor. /*! @@ -82,7 +82,17 @@ class BlowersMasel : public ArrheniusBase * @param w Average bond dissociation energy of the bond being formed and * broken in the reaction, in energy units [J/kmol] */ - BlowersMasel(double A, double b, double Ea0, double w); + BlowersMaselRate(double A, double b, double Ea0, double w); + + explicit BlowersMaselRate(const AnyMap& node, const UnitStack& rate_units={}) + : BlowersMaselRate() + { + setParameters(node, rate_units); + } + + unique_ptr newMultiRate() const override { + return unique_ptr(new MultiRate); + } virtual const std::string type() const override { return "Blowers-Masel"; @@ -177,8 +187,6 @@ class BlowersMasel : public ArrheniusBase double m_deltaH_R; //!< enthalpy change of reaction (in temperature units) }; -typedef BulkRate BlowersMaselRate; - } #endif diff --git a/include/cantera/kinetics/ChebyshevRate.h b/include/cantera/kinetics/ChebyshevRate.h index 37c66c4ad4c..21144eae08b 100644 --- a/include/cantera/kinetics/ChebyshevRate.h +++ b/include/cantera/kinetics/ChebyshevRate.h @@ -130,6 +130,8 @@ class ChebyshevRate final : public ReactionRate } void getParameters(AnyMap& rateNode) const; + virtual void validate(const std::string& equation, const Kinetics& kin); + //! Update information specific to reaction /*! * @param shared_data data shared by all reactions of a given type diff --git a/include/cantera/kinetics/Custom.h b/include/cantera/kinetics/Custom.h index 6b5192ddfbd..95127188d3f 100644 --- a/include/cantera/kinetics/Custom.h +++ b/include/cantera/kinetics/Custom.h @@ -52,6 +52,8 @@ class CustomFunc1Rate final : public ReactionRate void getParameters(AnyMap& rateNode, const Units& rate_units=Units(0.)) const; using ReactionRate::getParameters; + virtual void validate(const std::string& equation, const Kinetics& kin) override; + //! Update information specific to reaction /*! * @param shared_data data shared by all reactions of a given type diff --git a/include/cantera/kinetics/Falloff.h b/include/cantera/kinetics/Falloff.h index 37efc67ae3d..51ed84e1495 100644 --- a/include/cantera/kinetics/Falloff.h +++ b/include/cantera/kinetics/Falloff.h @@ -215,6 +215,7 @@ class FalloffRate : public ReactionRate } void check(const std::string& equation, const AnyMap& node); + virtual void validate(const std::string& equation, const Kinetics& kin); //! Get flag indicating whether negative A values are permitted bool allowNegativePreExponentialFactor() const { @@ -237,24 +238,24 @@ class FalloffRate : public ReactionRate } //! Get reaction rate in the low-pressure limit - Arrhenius3& lowRate() { + ArrheniusRate& lowRate() { return m_lowRate; } //! Set reaction rate in the low-pressure limit - void setLowRate(const Arrhenius3& low); + void setLowRate(const ArrheniusRate& low); //! Get reaction rate in the high-pressure limit - Arrhenius3& highRate() { + ArrheniusRate& highRate() { return m_highRate; } //! Set reaction rate in the high-pressure limit - void setHighRate(const Arrhenius3& high); + void setHighRate(const ArrheniusRate& high); protected: - Arrhenius3 m_lowRate; //!< The reaction rate in the low-pressure limit - Arrhenius3 m_highRate; //!< The reaction rate in the high-pressure limit + ArrheniusRate m_lowRate; //!< The reaction rate in the low-pressure limit + ArrheniusRate m_highRate; //!< The reaction rate in the high-pressure limit bool m_chemicallyActivated; //!< Flag labeling reaction as chemically activated bool m_negativeA_ok; //!< Flag indicating whether negative A values are permitted @@ -283,7 +284,7 @@ class LindemannRate final : public FalloffRate } LindemannRate( - const Arrhenius3& low, const Arrhenius3& high, const vector_fp& c) + const ArrheniusRate& low, const ArrheniusRate& high, const vector_fp& c) : LindemannRate() { m_lowRate = low; @@ -344,7 +345,7 @@ class TroeRate final : public FalloffRate setParameters(node, rate_units); } - TroeRate(const Arrhenius3& low, const Arrhenius3& high, const vector_fp& c) + TroeRate(const ArrheniusRate& low, const ArrheniusRate& high, const vector_fp& c) : TroeRate() { m_lowRate = low; @@ -446,7 +447,7 @@ class SriRate final : public FalloffRate setParameters(node, rate_units); } - SriRate(const Arrhenius3& low, const Arrhenius3& high, const vector_fp& c) + SriRate(const ArrheniusRate& low, const ArrheniusRate& high, const vector_fp& c) : SriRate() { m_lowRate = low; @@ -556,7 +557,7 @@ class TsangRate final : public FalloffRate setParameters(node, rate_units); } - TsangRate(const Arrhenius3& low, const Arrhenius3& high, const vector_fp& c) + TsangRate(const ArrheniusRate& low, const ArrheniusRate& high, const vector_fp& c) : TsangRate() { m_lowRate = low; diff --git a/include/cantera/kinetics/InterfaceRate.h b/include/cantera/kinetics/InterfaceRate.h index 7c21b328ed7..48eddd92dd9 100644 --- a/include/cantera/kinetics/InterfaceRate.h +++ b/include/cantera/kinetics/InterfaceRate.h @@ -86,7 +86,7 @@ struct InterfaceData : public BlowersMaselData * terms, where the parameters \f$ (a_k, E_k, m_k) \f$ describe the dependency on the * surface coverage of species \f$ k, \theta_k \f$. The InterfaceRateBase class implements * terms related to coverage only, which allows for combinations with arbitrary rate - * parameterizations (for example Arrhenius and BlowersMasel). + * parameterizations (for example Arrhenius and BlowersMaselRate). */ class InterfaceRateBase { @@ -359,9 +359,12 @@ class InterfaceRate : public RateType, public InterfaceRateBase using RateType::RateType; // inherit constructors //! Constructor based on AnyMap content - InterfaceRate(const AnyMap& node, const UnitStack& rate_units={}) { + InterfaceRate(const AnyMap& node, const UnitStack& rate_units) { setParameters(node, rate_units); } + explicit InterfaceRate(const AnyMap& node) { + setParameters(node, {}); + } unique_ptr newMultiRate() const override { return unique_ptr( @@ -377,25 +380,12 @@ class InterfaceRate : public RateType, public InterfaceRateBase const AnyMap& node, const UnitStack& rate_units) override { InterfaceRateBase::setParameters(node); - RateType::m_negativeA_ok = node.getBool("negative-A", false); - if (!node.hasKey("rate-constant")) { - RateType::setRateParameters(AnyValue(), node.units(), rate_units); - return; - } - RateType::setRateParameters(node["rate-constant"], node.units(), rate_units); + RateType::setParameters(node, rate_units); } virtual void getParameters(AnyMap& node) const override { + RateType::getParameters(node); node["type"] = type(); - if (RateType::m_negativeA_ok) { - node["negative-A"] = true; - } - AnyMap rateNode; - RateType::getRateParameters(rateNode); - if (!rateNode.empty()) { - // RateType object is configured - node["rate-constant"] = std::move(rateNode); - } InterfaceRateBase::getParameters(node); } @@ -455,8 +445,8 @@ class InterfaceRate : public RateType, public InterfaceRateBase } }; -using InterfaceArrheniusRate = InterfaceRate; -using InterfaceBlowersMaselRate = InterfaceRate; +using InterfaceArrheniusRate = InterfaceRate; +using InterfaceBlowersMaselRate = InterfaceRate; //! A class template for interface sticking rate specifications @@ -470,7 +460,11 @@ class StickingRate : public RateType, public StickingCoverage using RateType::RateType; // inherit constructors //! Constructor based on AnyMap content - StickingRate(const AnyMap& node, const UnitStack& rate_units={}) { + StickingRate(const AnyMap& node, const UnitStack& rate_units) { + // sticking coefficients are dimensionless + setParameters(node, Units(1.0)); + } + explicit StickingRate(const AnyMap& node) { // sticking coefficients are dimensionless setParameters(node, Units(1.0)); } @@ -521,7 +515,7 @@ class StickingRate : public RateType, public StickingCoverage } virtual void validate(const std::string &equation, const Kinetics& kin) override { - ReactionRate::validate(equation, kin); + RateType::validate(equation, kin); fmt::memory_buffer err_reactions; double T[] = {200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0}; for (size_t i=0; i < 6; i++) { @@ -594,8 +588,8 @@ class StickingRate : public RateType, public StickingCoverage } }; -using StickingArrheniusRate = StickingRate; -using StickingBlowersMaselRate = StickingRate; +using StickingArrheniusRate = StickingRate; +using StickingBlowersMaselRate = StickingRate; } #endif diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 44da8ef3e5a..29fb438ec83 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -621,7 +621,7 @@ class Kinetics /** * Calculate derivatives for forward rate constants with respect to temperature * at constant pressure, molar concentration and mole fractions - * @param dkfwd[out] Output vector of derivatives. Length: nReactions(). + * @param[out] dkfwd Output vector of derivatives. Length: nReactions(). */ virtual void getFwdRateConstants_ddT(double* dkfwd) { @@ -632,7 +632,7 @@ class Kinetics /** * Calculate derivatives for forward rate constants with respect to pressure * at constant temperature, molar concentration and mole fractions. - * @param dkfwd[out] Output vector of derivatives. Length: nReactions(). + * @param[out] dkfwd Output vector of derivatives. Length: nReactions(). */ virtual void getFwdRateConstants_ddP(double* dkfwd) { @@ -643,7 +643,7 @@ class Kinetics /** * Calculate derivatives for forward rate constants with respect to molar * concentration at constant temperature, pressure and mole fractions. - * @param dkfwd[out] Output vector of derivatives. Length: nReactions(). + * @param[out] dkfwd Output vector of derivatives. Length: nReactions(). * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -657,7 +657,7 @@ class Kinetics /** * Calculate derivatives for forward rates-of-progress with respect to temperature * at constant pressure, molar concentration and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). */ virtual void getFwdRatesOfProgress_ddT(double* drop) { @@ -668,7 +668,7 @@ class Kinetics /** * Calculate derivatives for forward rates-of-progress with respect to pressure * at constant temperature, molar concentration and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). */ virtual void getFwdRatesOfProgress_ddP(double* drop) { @@ -679,7 +679,7 @@ class Kinetics /** * Calculate derivatives for forward rates-of-progress with respect to molar * concentration at constant temperature, pressure and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -695,8 +695,8 @@ class Kinetics * mole fractions at constant temperature, pressure and molar concentration. * * The method returns a matrix with nReactions rows and nTotalSpecies columns. - * For a derivative with respect to \f$X_i$, all other \f$X_j$ are held constant, - * rather than enforcing \f$\sum X_j = 1$. + * For a derivative with respect to \f$X_i\f$, all other \f$X_j\f$ are held + * constant, rather than enforcing \f$\sum X_j = 1\f$. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -710,7 +710,7 @@ class Kinetics /** * Calculate derivatives for reverse rates-of-progress with respect to temperature * at constant pressure, molar concentration and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). */ virtual void getRevRatesOfProgress_ddT(double* drop) { @@ -721,7 +721,7 @@ class Kinetics /** * Calculate derivatives for reverse rates-of-progress with respect to pressure * at constant temperature, molar concentration and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). */ virtual void getRevRatesOfProgress_ddP(double* drop) { @@ -732,7 +732,7 @@ class Kinetics /** * Calculate derivatives for reverse rates-of-progress with respect to molar * concentration at constant temperature, pressure and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -748,8 +748,8 @@ class Kinetics * mole fractions at constant temperature, pressure and molar concentration. * * The method returns a matrix with nReactions rows and nTotalSpecies columns. - * For a derivative with respect to \f$X_i$, all other \f$X_j$ are held constant, - * rather than enforcing \f$\sum X_j = 1$. + * For a derivative with respect to \f$X_i\f$, all other \f$X_j\f$ are held + * constant, rather than enforcing \f$\sum X_j = 1\f$. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -763,7 +763,7 @@ class Kinetics /** * Calculate derivatives for net rates-of-progress with respect to temperature * at constant pressure, molar concentration and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). */ virtual void getNetRatesOfProgress_ddT(double* drop) { @@ -774,7 +774,7 @@ class Kinetics /** * Calculate derivatives for net rates-of-progress with respect to pressure * at constant temperature, molar concentration and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). */ virtual void getNetRatesOfProgress_ddP(double* drop) { @@ -785,7 +785,7 @@ class Kinetics /** * Calculate derivatives for net rates-of-progress with respect to molar * concentration at constant temperature, pressure and mole fractions. - * @param drop[out] Output vector of derivatives. Length: nReactions(). + * @param[out] drop Output vector of derivatives. Length: nReactions(). * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -801,8 +801,8 @@ class Kinetics * mole fractions at constant temperature, pressure and molar concentration. * * The method returns a matrix with nReactions rows and nTotalSpecies columns. - * For a derivative with respect to \f$X_i$, all other \f$X_j$ are held constant, - * rather than enforcing \f$\sum X_j = 1$. + * For a derivative with respect to \f$X_i\f$, all other \f$X_j\f$ are held + * constant, rather than enforcing \f$\sum X_j = 1\f$. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -816,21 +816,21 @@ class Kinetics /** * Calculate derivatives for species creation rates with respect to temperature * at constant pressure, molar concentration and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. */ void getCreationRates_ddT(double* dwdot); /** * Calculate derivatives for species creation rates with respect to pressure * at constant temperature, molar concentration and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. */ void getCreationRates_ddP(double* dwdot); /** * Calculate derivatives for species creation rates with respect to molar * concentration at constant temperature, pressure and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -842,8 +842,8 @@ class Kinetics * mole fractions at constant temperature, pressure and molar concentration. * * The method returns a matrix with nTotalSpecies rows and nTotalSpecies columns. - * For a derivative with respect to \f$X_i$, all other \f$X_j$ are held constant, - * rather than enforcing \f$\sum X_j = 1$. + * For a derivative with respect to \f$X_i\f$, all other \f$X_j\f$ are held + * constant, rather than enforcing \f$\sum X_j = 1\f$. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -853,21 +853,21 @@ class Kinetics /** * Calculate derivatives for species destruction rates with respect to temperature * at constant pressure, molar concentration and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. */ void getDestructionRates_ddT(double* dwdot); /** * Calculate derivatives for species destruction rates with respect to pressure * at constant temperature, molar concentration and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. */ void getDestructionRates_ddP(double* dwdot); /** * Calculate derivatives for species destruction rates with respect to molar * concentration at constant temperature, pressure and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -879,8 +879,8 @@ class Kinetics * mole fractions at constant temperature, pressure and molar concentration. * * The method returns a matrix with nTotalSpecies rows and nTotalSpecies columns. - * For a derivative with respect to \f$X_i$, all other \f$X_j$ are held constant, - * rather than enforcing \f$\sum X_j = 1$. + * For a derivative with respect to \f$X_i\f$, all other \f$X_j\f$ are held + * constant, rather than enforcing \f$\sum X_j = 1\f$. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -890,21 +890,21 @@ class Kinetics /** * Calculate derivatives for species net production rates with respect to * temperature at constant pressure, molar concentration and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. */ void getNetProductionRates_ddT(double* dwdot); /** * Calculate derivatives for species net production rates with respect to pressure * at constant temperature, molar concentration and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. */ void getNetProductionRates_ddP(double* dwdot); /** * Calculate derivatives for species net production rates with respect to molar * concentration at constant temperature, pressure and mole fractions. - * @param dwdot[out] Output vector of derivatives. Length: m_kk. + * @param[out] dwdot Output vector of derivatives. Length: m_kk. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. @@ -916,8 +916,8 @@ class Kinetics * mole fractions at constant temperature, pressure and molar concentration. * * The method returns a matrix with nTotalSpecies rows and nTotalSpecies columns. - * For a derivative with respect to \f$X_i$, all other \f$X_j$ are held constant, - * rather than enforcing \f$\sum X_j = 1$. + * For a derivative with respect to \f$X_i\f$, all other \f$X_j\f$ are held constant, + * rather than enforcing \f$\sum X_j = 1\f$. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. diff --git a/include/cantera/kinetics/PlogRate.h b/include/cantera/kinetics/PlogRate.h index 128a979f390..f9564db13f0 100644 --- a/include/cantera/kinetics/PlogRate.h +++ b/include/cantera/kinetics/PlogRate.h @@ -81,7 +81,7 @@ class PlogRate final : public ReactionRate PlogRate(); //! Constructor from Arrhenius rate expressions at a set of pressures - explicit PlogRate(const std::multimap& rates); + explicit PlogRate(const std::multimap& rates); //! Constructor using legacy Arrhenius2 framework explicit PlogRate(const std::multimap& rates); @@ -104,15 +104,6 @@ class PlogRate final : public ReactionRate */ void setParameters(const AnyMap& node, const UnitStack& units); - //! Perform object setup based on reaction rate information - /*! - * @param rates vector of AnyMap containing rate information - * @param units unit system - * @param rate_units unit definitions specific to rate information - */ - void setRateParameters(const std::vector& rates, - const UnitSystem& units, const Units& rate_units); - void getParameters(AnyMap& rateNode, const Units& rate_units) const; void getParameters(AnyMap& rateNode) const { return getParameters(rateNode, Units(0)); @@ -143,7 +134,7 @@ class PlogRate final : public ReactionRate void setup(const std::multimap& rates); //! Set up Plog object - void setRates(const std::multimap& rates); + void setRates(const std::multimap& rates); //! Update concentration-dependent parts of the rate coefficient. //! @param c natural log of the pressure in Pa @@ -226,14 +217,14 @@ class PlogRate final : public ReactionRate //! Return the pressures and Arrhenius expressions which comprise this //! reaction. - std::multimap getRates() const; + std::multimap getRates() const; protected: //! log(p) to (index range) in the rates_ vector std::map> pressures_; // Rate expressions which are referenced by the indices stored in pressures_ - std::vector rates_; + std::vector rates_; double logP_; //!< log(p) at the current state double logP1_, logP2_; //!< log(p) at the lower / upper pressure reference diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 816ca97bb77..242a62bc854 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -36,9 +36,9 @@ class Func1; * k_f = A T^b \exp (-E/RT) * \f] * - * @deprecated To be removed after Cantera 2.6. See Arrhenius3 / ArrheniusRate. + * @deprecated To be removed after Cantera 2.6. See ArrheniusRate. */ -class Arrhenius2 final : public Arrhenius3 +class Arrhenius2 final : public ArrheniusRate { public: //! Default constructor. @@ -60,15 +60,15 @@ class Arrhenius2 final : public Arrhenius3 const UnitSystem& units, const Units& rate_units); //! Converting constructor (to facilitate back-ward compatibility) - Arrhenius2(const Arrhenius3& other); + Arrhenius2(const ArrheniusRate& other); void setRateParameters(const AnyValue& rate, const UnitSystem& units, const Units& rate_units); - using Arrhenius3::setRateParameters; + using ArrheniusRate::setRateParameters; //! Return parameters - two-parameter version void getParameters(AnyMap& node, const Units& rate_units) const; - using Arrhenius3::getParameters; + using ArrheniusRate::getParameters; //! Update concentration-dependent parts of the rate coefficient. /*! @@ -202,7 +202,7 @@ class SurfaceArrhenius #ifdef CT_NO_LEGACY_REACTIONS_26 -typedef Arrhenius3 Arrhenius; +typedef ArrheniusRate Arrhenius; #else typedef Arrhenius2 Arrhenius; #endif diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 57cb567e41a..25f96d93289 100644 --- a/include/cantera/kinetics/TwoTempPlasmaRate.h +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -58,10 +58,10 @@ struct TwoTempPlasmaData : public ReactionData * * @ingroup arrheniusGroup */ -class TwoTempPlasma : public ArrheniusBase +class TwoTempPlasmaRate : public ArrheniusBase { public: - TwoTempPlasma(); + TwoTempPlasmaRate(); //! Constructor. /*! @@ -71,7 +71,16 @@ class TwoTempPlasma : public ArrheniusBase * @param Ea Activation energy in energy units [J/kmol] * @param EE Activation electron energy in energy units [J/kmol] */ - TwoTempPlasma(double A, double b, double Ea=0.0, double EE=0.0); + TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); + + TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}) : TwoTempPlasmaRate() { + setParameters(node, rate_units); + } + + unique_ptr newMultiRate() const override { + return unique_ptr( + new MultiRate); + } virtual const std::string type() const override { return "two-temperature-plasma"; @@ -106,8 +115,6 @@ class TwoTempPlasma : public ArrheniusBase } }; -typedef BulkRate TwoTempPlasmaRate; - } #endif diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index e51aafe36ec..16e0b4041a8 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -462,17 +462,14 @@ cdef extern from "cantera/kinetics/Arrhenius.h" namespace "Cantera": cbool allowNegativePreExponentialFactor() void setAllowNegativePreExponentialFactor(bool) - cdef cppclass CxxArrhenius "Cantera::Arrhenius3" (CxxArrheniusBase): - CxxArrhenius(double, double, double) + cdef cppclass CxxArrheniusRate "Cantera::ArrheniusRate" (CxxArrheniusBase): + CxxArrheniusRate(double, double, double) + CxxArrheniusRate(CxxAnyMap) except +translate_exception double evalRate(double, double) - cdef cppclass CxxArrhenius2 "Cantera::Arrhenius2" (CxxArrhenius): + cdef cppclass CxxArrhenius2 "Cantera::Arrhenius2" (CxxArrheniusRate): CxxArrhenius2(double, double, double) - cdef cppclass CxxArrheniusRate "Cantera::ArrheniusRate" (CxxArrhenius): - CxxArrheniusRate(CxxAnyMap) except +translate_exception - CxxArrheniusRate(double, double, double) - cdef cppclass CxxBlowersMasel "Cantera::BlowersMasel" (CxxArrheniusBase): CxxBlowersMasel(double, double, double, double) double evalRate(double, double) @@ -504,10 +501,10 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": void setAllowNegativePreExponentialFactor(bool) cbool chemicallyActivated() void setChemicallyActivated(bool) - CxxArrhenius& lowRate() - void setLowRate(CxxArrhenius&) except +translate_exception - CxxArrhenius& highRate() - void setHighRate(CxxArrhenius&) except +translate_exception + CxxArrheniusRate& lowRate() + void setLowRate(CxxArrheniusRate&) except +translate_exception + CxxArrheniusRate& highRate() + void setHighRate(CxxArrheniusRate&) except +translate_exception void getFalloffCoeffs(vector[double]&) void setFalloffCoeffs(vector[double]&) except +translate_exception double evalF(double, double) except +translate_exception @@ -527,8 +524,8 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxPlogRate "Cantera::PlogRate" (CxxReactionRate): CxxPlogRate() CxxPlogRate(CxxAnyMap) except +translate_exception - CxxPlogRate(multimap[double, CxxArrhenius]) - multimap[double, CxxArrhenius] getRates() + CxxPlogRate(multimap[double, CxxArrheniusRate]) + multimap[double, CxxArrheniusRate] getRates() cdef cppclass CxxChebyshevRate "Cantera::ChebyshevRate" (CxxReactionRate): CxxChebyshevRate() @@ -565,12 +562,12 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": double stickingWeight() void setStickingWeight(double) - cdef cppclass CxxInterfaceArrheniusRate "Cantera::InterfaceArrheniusRate" (CxxReactionRate, CxxArrhenius, CxxInterfaceRateBase): + cdef cppclass CxxInterfaceArrheniusRate "Cantera::InterfaceArrheniusRate" (CxxReactionRate, CxxArrheniusRate, CxxInterfaceRateBase): CxxInterfaceArrheniusRate() CxxInterfaceArrheniusRate(CxxAnyMap) except +translate_exception CxxInterfaceArrheniusRate(double, double, double) - cdef cppclass CxxStickingArrheniusRate "Cantera::StickingArrheniusRate" (CxxReactionRate, CxxArrhenius, CxxStickingCoverage): + cdef cppclass CxxStickingArrheniusRate "Cantera::StickingArrheniusRate" (CxxReactionRate, CxxArrheniusRate, CxxStickingCoverage): CxxStickingArrheniusRate() CxxStickingArrheniusRate(CxxAnyMap) except +translate_exception CxxStickingArrheniusRate(double, double, double) @@ -1426,11 +1423,11 @@ cdef class CustomReaction(Reaction): cdef class Arrhenius: cdef CxxArrhenius2* legacy # used by legacy objects only - cdef CxxArrhenius* base + cdef CxxArrheniusRate* base cdef cbool own_rate cdef Reaction reaction # parent reaction, to prevent garbage collection @staticmethod - cdef wrap(CxxArrhenius*) + cdef wrap(CxxArrheniusRate*) cdef class Falloff: cdef shared_ptr[CxxFalloff] _falloff diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index 8737fd418dd..11507b7f3c2 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -216,7 +216,7 @@ cdef class ArrheniusRate(ArrheniusRateBase): cdef set_cxx_object(self): self.rate = self._rate.get() - self.base = self.rate + self.base = self.rate cdef CxxArrheniusRate* cxx_object(self): return self.rate @@ -528,17 +528,17 @@ cdef class PlogRate(ReactionRate): """ def __get__(self): rates = [] - cdef multimap[double, CxxArrhenius] cxxrates - cdef pair[double, CxxArrhenius] p_rate + cdef multimap[double, CxxArrheniusRate] cxxrates + cdef pair[double, CxxArrheniusRate] p_rate cxxrates = self.cxx_object().getRates() for p_rate in cxxrates: rates.append((p_rate.first, copyArrhenius(&p_rate.second))) return rates def __set__(self, rates): - cdef multimap[double, CxxArrhenius] ratemap + cdef multimap[double, CxxArrheniusRate] ratemap cdef Arrhenius rate - cdef pair[double, CxxArrhenius] item + cdef pair[double, CxxArrheniusRate] item for p, rate in rates: item.first = p item.second = deref(rate.base) @@ -782,7 +782,7 @@ cdef class InterfaceArrheniusRate(InterfaceRateBase): cdef set_cxx_object(self): self.rate = self._rate.get() - self.base = self.rate + self.base = self.rate self.interface = self.cxx_object() cdef CxxInterfaceArrheniusRate* cxx_object(self): @@ -810,7 +810,7 @@ cdef class InterfaceBlowersMaselRate(InterfaceRateBase): cdef set_cxx_object(self): self.rate = self._rate.get() - self.base = self.rate + self.base = self.rate self.interface = self.cxx_object() cdef CxxInterfaceBlowersMaselRate* cxx_object(self): @@ -920,7 +920,7 @@ cdef class StickingArrheniusRate(StickRateBase): cdef set_cxx_object(self): self.rate = self._rate.get() - self.base = self.rate + self.base = self.rate self.stick = self.cxx_object() self.interface = self.stick @@ -947,7 +947,7 @@ cdef class StickingBlowersMaselRate(StickRateBase): cdef set_cxx_object(self): self.rate = self._rate.get() - self.base = self.rate + self.base = self.rate self.stick = self.cxx_object() self.interface = self.stick @@ -1920,7 +1920,7 @@ cdef class Arrhenius: """ def __cinit__(self, A=0, b=0, E=0, init=True): if init: - self.base = new CxxArrhenius(A, b, E) + self.base = new CxxArrheniusRate(A, b, E) self.own_rate = True self.reaction = None else: @@ -1931,7 +1931,7 @@ cdef class Arrhenius: del self.base @staticmethod - cdef wrap(CxxArrhenius* rate): + cdef wrap(CxxArrheniusRate* rate): r = Arrhenius(init=False) r.base = rate r.reaction = None @@ -1968,7 +1968,7 @@ cdef class Arrhenius: return self.base.evalRate(np.log(T), 1/T) -cdef wrapArrhenius(CxxArrhenius* rate, Reaction reaction): +cdef wrapArrhenius(CxxArrheniusRate* rate, Reaction reaction): r = Arrhenius(init=False) r.base = rate if reaction.uses_legacy: @@ -1982,7 +1982,7 @@ cdef copyArrhenius2(CxxArrhenius2* rate): rate.activationEnergy()) return r -cdef copyArrhenius(CxxArrhenius* rate): +cdef copyArrhenius(CxxArrheniusRate* rate): r = Arrhenius(rate.preExponentialFactor(), rate.temperatureExponent(), rate.activationEnergy()) return r diff --git a/interfaces/cython/cantera/test/test_reaction.py b/interfaces/cython/cantera/test/test_reaction.py index b198f9ae617..62d0ccf5cfa 100644 --- a/interfaces/cython/cantera/test/test_reaction.py +++ b/interfaces/cython/cantera/test/test_reaction.py @@ -1119,18 +1119,19 @@ def test_no_rate(self): else: self.assertIsNaN(self.eval_rate(rxn.rate)) - sol2 = ct.Solution(thermo=self.soln.thermo_model, kinetics=self.soln.kinetics_model, - species=self.species, reactions=[rxn], adjacent=self.adj) - sol2.TPX = self.soln.TPX if self._legacy: + sol2 = ct.Solution(thermo=self.soln.thermo_model, + kinetics=self.soln.kinetics_model, + species=self.species, + reactions=[rxn], adjacent=self.adj) + sol2.TPX = self.soln.TPX self.assertNear(sol2.forward_rate_constants[0], 0.) self.assertNear(sol2.net_rates_of_progress[0], 0.) - elif not ct.debug_mode_enabled(): - self.assertIsNaN(sol2.forward_rate_constants[0]) - self.assertIsNaN(sol2.net_rates_of_progress[0]) else: - with self.assertRaisesRegex(ct.CanteraError, "not finite"): - sol2.net_rates_of_progress + with self.assertRaisesRegex(ct.CanteraError, "validate"): + ct.Solution(thermo=self.soln.thermo_model, + kinetics=self.soln.kinetics_model, + species=self.species, reactions=[rxn], adjacent=self.adj) def test_replace_rate(self): # check replacing reaction rate expression diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index ebb6c55de4d..483a5f084aa 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -108,7 +108,32 @@ void ArrheniusBase::getRateParameters(AnyMap& node) const node[m_E4_str].setQuantity(m_E4_R, "K", true); } node.setFlowStyle(); +} +void ArrheniusBase::setParameters(const AnyMap& node, const UnitStack& rate_units) +{ + ReactionRate::setParameters(node, rate_units); + m_negativeA_ok = node.getBool("negative-A", false); + if (!node.hasKey("rate-constant")) { + setRateParameters(AnyValue(), node.units(), rate_units); + return; + } + setRateParameters(node["rate-constant"], node.units(), rate_units); +} + +void ArrheniusBase::getParameters(AnyMap& node) const { + if (m_negativeA_ok) { + node["negative-A"] = true; + } + AnyMap rateNode; + getRateParameters(rateNode); + if (!rateNode.empty()) { + // RateType object is configured + node["rate-constant"] = std::move(rateNode); + } + if (type() != "Arrhenius") { + node["type"] = type(); + } } void ArrheniusBase::check(const std::string& equation, const AnyMap& node) @@ -126,6 +151,14 @@ void ArrheniusBase::check(const std::string& equation, const AnyMap& node) } } +void ArrheniusBase::validate(const std::string& equation, const Kinetics& kin) +{ + if (isnan(m_A) || isnan(m_b)) { + throw CanteraError("ArrheniusBase::validate", + "Rate object for reaction '{}' is not configured.", equation); + } +} + bool ArrheniusData::update(const ThermoPhase& phase, const Kinetics& kin) { double T = phase.temperature(); diff --git a/src/kinetics/BlowersMaselRate.cpp b/src/kinetics/BlowersMaselRate.cpp index 81ed5f2211e..8cfa5738acd 100644 --- a/src/kinetics/BlowersMaselRate.cpp +++ b/src/kinetics/BlowersMaselRate.cpp @@ -43,14 +43,14 @@ bool BlowersMaselData::update(const ThermoPhase& phase, const Kinetics& kin) return changed; } -BlowersMasel::BlowersMasel() +BlowersMaselRate::BlowersMaselRate() : m_deltaH_R(0.) { m_Ea_str = "Ea0"; m_E4_str = "w"; } -BlowersMasel::BlowersMasel(double A, double b, double Ea0, double w) +BlowersMaselRate::BlowersMaselRate(double A, double b, double Ea0, double w) : ArrheniusBase(A, b, Ea0) , m_deltaH_R(0.) { @@ -59,15 +59,15 @@ BlowersMasel::BlowersMasel(double A, double b, double Ea0, double w) m_E4_R = w / GasConstant; } -double BlowersMasel::ddTScaledFromStruct(const BlowersMaselData& shared_data) const +double BlowersMaselRate::ddTScaledFromStruct(const BlowersMaselData& shared_data) const { - warn_user("BlowersMasel::ddTScaledFromStruct", + warn_user("BlowersMaselRate::ddTScaledFromStruct", "Temperature derivative does not consider changes of reaction enthalpy."); double Ea_R = effectiveActivationEnergy_R(m_deltaH_R); return (Ea_R * shared_data.recipT + m_b) * shared_data.recipT; } -void BlowersMasel::setContext(const Reaction& rxn, const Kinetics& kin) +void BlowersMaselRate::setContext(const Reaction& rxn, const Kinetics& kin) { m_stoich_coeffs.clear(); for (const auto& sp : rxn.reactants) { diff --git a/src/kinetics/ChebyshevRate.cpp b/src/kinetics/ChebyshevRate.cpp index dc7b94bd17d..61fbcc06eb1 100644 --- a/src/kinetics/ChebyshevRate.cpp +++ b/src/kinetics/ChebyshevRate.cpp @@ -171,4 +171,12 @@ void ChebyshevRate::getParameters(AnyMap& rateNode) const rateNode["data"].setQuantity(coeffs, converter); } +void ChebyshevRate::validate(const std::string& equation, const Kinetics& kin) +{ + if (m_coeffs.data().empty() || isnan(m_coeffs(0, 0))) { + throw CanteraError("ChebyshevRate::validate", + "Rate object for reaction '{}' is not configured.", equation); + } +} + } diff --git a/src/kinetics/Custom.cpp b/src/kinetics/Custom.cpp index 52ffc684dce..fc935f3ebc5 100644 --- a/src/kinetics/Custom.cpp +++ b/src/kinetics/Custom.cpp @@ -19,6 +19,14 @@ void CustomFunc1Rate::setRateFunction(shared_ptr f) m_ratefunc = f; } +void CustomFunc1Rate::validate(const std::string& equation, const Kinetics& kin) +{ + if (!m_ratefunc) { + throw CanteraError("CustomFunc1Rate::validate", + "Rate object for reaction '{}' is not configured.", equation); + } +} + double CustomFunc1Rate::evalFromStruct(const ArrheniusData& shared_data) const { if (m_ratefunc) { diff --git a/src/kinetics/Falloff.cpp b/src/kinetics/Falloff.cpp index 414da9232b6..d383f4547ca 100644 --- a/src/kinetics/Falloff.cpp +++ b/src/kinetics/Falloff.cpp @@ -87,9 +87,9 @@ void FalloffRate::init(const vector_fp& c) setFalloffCoeffs(c); } -void FalloffRate::setLowRate(const Arrhenius3& low) +void FalloffRate::setLowRate(const ArrheniusRate& low) { - Arrhenius3 _low = low; + ArrheniusRate _low = low; _low.setAllowNegativePreExponentialFactor(m_negativeA_ok); _low.check("", AnyMap()); if (_low.preExponentialFactor() * m_highRate.preExponentialFactor() < 0.) { @@ -100,9 +100,9 @@ void FalloffRate::setLowRate(const Arrhenius3& low) m_lowRate = std::move(_low); } -void FalloffRate::setHighRate(const Arrhenius3& high) +void FalloffRate::setHighRate(const ArrheniusRate& high) { - Arrhenius3 _high = high; + ArrheniusRate _high = high; _high.setAllowNegativePreExponentialFactor(m_negativeA_ok); _high.check("", AnyMap()); if (m_lowRate.preExponentialFactor() * _high.preExponentialFactor() < 0.) { @@ -149,12 +149,12 @@ void FalloffRate::setParameters(const AnyMap& node, const UnitStack& rate_units) } } if (node.hasKey("low-P-rate-constant")) { - m_lowRate = Arrhenius3( + m_lowRate = ArrheniusRate( node["low-P-rate-constant"], node.units(), low_rate_units); m_lowRate.setAllowNegativePreExponentialFactor(m_negativeA_ok); } if (node.hasKey("high-P-rate-constant")) { - m_highRate = Arrhenius3( + m_highRate = ArrheniusRate( node["high-P-rate-constant"], node.units(), high_rate_units); m_highRate.setAllowNegativePreExponentialFactor(m_negativeA_ok); } @@ -200,6 +200,12 @@ void FalloffRate::check(const std::string& equation, const AnyMap& node) } } +void FalloffRate::validate(const std::string& equation, const Kinetics& kin) +{ + m_lowRate.validate(equation, kin); + m_highRate.validate(equation, kin); +} + void TroeRate::setFalloffCoeffs(const vector_fp& c) { if (c.size() != 3 && c.size() != 4) { @@ -314,10 +320,6 @@ void TroeRate::getParameters(AnyMap& node) const if (std::abs(m_t2) > SmallNumber) { params["T2"] = m_t2; } - // This can't be converted to a different unit system because the dimensions of - // the rate constant were not set. Can occur if the reaction was created outside - // the context of a Kinetics object and never added to a Kinetics object. - node["__unconvertible__"] = true; } params.setFlowStyle(); node["Troe"] = std::move(params); @@ -434,10 +436,6 @@ void SriRate::getParameters(AnyMap& node) const params["D"] = m_d; params["E"] = m_e; } - // This can't be converted to a different unit system because the dimensions of - // the rate constant were not set. Can occur if the reaction was created outside - // the context of a Kinetics object and never added to a Kinetics object. - node["__unconvertible__"] = true; } params.setFlowStyle(); node["SRI"] = std::move(params); diff --git a/src/kinetics/InterfaceRate.cpp b/src/kinetics/InterfaceRate.cpp index 842c8dc8c40..b6a4c103c36 100644 --- a/src/kinetics/InterfaceRate.cpp +++ b/src/kinetics/InterfaceRate.cpp @@ -1,4 +1,4 @@ -//! @file Coverage.cpp +//! @file InterfaceRate.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. diff --git a/src/kinetics/PlogRate.cpp b/src/kinetics/PlogRate.cpp index 44c5a51e566..397b654e3c0 100644 --- a/src/kinetics/PlogRate.cpp +++ b/src/kinetics/PlogRate.cpp @@ -56,7 +56,7 @@ PlogRate::PlogRate() { } -PlogRate::PlogRate(const std::multimap& rates) +PlogRate::PlogRate(const std::multimap& rates) : PlogRate() { setRates(rates); @@ -70,29 +70,17 @@ PlogRate::PlogRate(const std::multimap& rates) void PlogRate::setParameters(const AnyMap& node, const UnitStack& units) { - auto rate_units = units.product(); + std::multimap multi_rates; if (!node.hasKey("rate-constants")) { - PlogRate::setRateParameters(std::vector (), node.units(), rate_units); - return; - } - - setRateParameters(node.at("rate-constants").asVector(), - node.units(), rate_units); -} - -void PlogRate::setRateParameters(const std::vector& rates, - const UnitSystem& units, const Units& rate_units) -{ - std::multimap multi_rates; - if (rates.size()) { + // ensure that reaction rate can be evaluated (but returns NaN) + multi_rates.insert({1.e-7, ArrheniusRate(NAN, NAN, NAN)}); + multi_rates.insert({1.e14, ArrheniusRate(NAN, NAN, NAN)}); + } else { + auto& rates = node["rate-constants"].asVector(); for (const auto& rate : rates) { multi_rates.insert({rate.convert("P", "Pa"), - Arrhenius3(AnyValue(rate), units, rate_units)}); + ArrheniusRate(AnyValue(rate), node.units(), units)}); } - } else { - // ensure that reaction rate can be evaluated (but returns NaN) - multi_rates.insert({1.e-7, Arrhenius3(NAN, NAN, NAN)}); - multi_rates.insert({1.e14, Arrhenius3(NAN, NAN, NAN)}); } setRates(multi_rates); } @@ -123,14 +111,14 @@ void PlogRate::getParameters(AnyMap& rateNode, const Units& rate_units) const void PlogRate::setup(const std::multimap& rates) { - std::multimap rates2; + std::multimap rates2; for (const auto& item : rates) { rates2.emplace(item.first, item.second); } setRates(rates2); } -void PlogRate::setRates(const std::multimap& rates) +void PlogRate::setRates(const std::multimap& rates) { size_t j = 0; rates_.clear(); @@ -168,7 +156,7 @@ void PlogRate::validate(const std::string& equation) for (size_t p = ilow1_; p < ilow2_; p++) { k += rates_.at(p).evalRate(log(T[i]), 1.0 / T[i]); } - if (k < 0) { + if (!(k > 0)) { fmt_append(err_reactions, "\nInvalid rate coefficient for reaction '{}'\n" "at P = {:.5g}, T = {:.1f}\n", @@ -191,9 +179,9 @@ std::vector > PlogRate::rates() const return out; } -std::multimap PlogRate::getRates() const +std::multimap PlogRate::getRates() const { - std::multimap rateMap; + std::multimap rateMap; // initial preincrement to skip rate for P --> 0 for (auto iter = ++pressures_.begin(); iter->first < 1000; // skip rates for (P --> infinity) diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 151982dfb99..6266b62ff39 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -9,7 +9,7 @@ namespace Cantera { Arrhenius2::Arrhenius2() - : Arrhenius3() + : ArrheniusRate() { m_b = 0.0; m_A = 0.0; @@ -17,7 +17,7 @@ Arrhenius2::Arrhenius2() } Arrhenius2::Arrhenius2(doublereal A, doublereal b, doublereal E) - : Arrhenius3(A, b, E * GasConstant) + : ArrheniusRate(A, b, E * GasConstant) { if (m_A <= 0.0) { m_logA = -1.0E300; @@ -30,8 +30,8 @@ Arrhenius2::Arrhenius2(const AnyValue& rate, setRateParameters(rate, units, rate_units); } -Arrhenius2::Arrhenius2(const Arrhenius3& other) - : Arrhenius3(other.preExponentialFactor(), +Arrhenius2::Arrhenius2(const ArrheniusRate& other) + : ArrheniusRate(other.preExponentialFactor(), other.temperatureExponent(), other.activationEnergy()) { @@ -41,7 +41,7 @@ void Arrhenius2::setRateParameters(const AnyValue& rate, const UnitSystem& units, const Units& rate_units) { UnitStack units_stack(rate_units); - Arrhenius3::setRateParameters(rate, units, units_stack); + ArrheniusRate::setRateParameters(rate, units, units_stack); if (m_A <= 0.0) { m_logA = -1.0E300; } diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index 7b5a2a0f373..4b75888abff 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -45,14 +45,14 @@ void TwoTempPlasmaData::updateTe(double Te) recipTe = 1./Te; } -TwoTempPlasma::TwoTempPlasma() +TwoTempPlasmaRate::TwoTempPlasmaRate() : ArrheniusBase() { m_Ea_str = "Ea-gas"; m_E4_str = "Ea-electron"; } -TwoTempPlasma::TwoTempPlasma(double A, double b, double Ea, double EE) +TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) : ArrheniusBase(A, b, Ea) { m_Ea_str = "Ea-gas"; @@ -60,20 +60,20 @@ TwoTempPlasma::TwoTempPlasma(double A, double b, double Ea, double EE) m_E4_R = EE / GasConstant; } -double TwoTempPlasma::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const +double TwoTempPlasmaRate::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const { - warn_user("TwoTempPlasma::ddTScaledFromStruct", + warn_user("TwoTempPlasmaRate::ddTScaledFromStruct", "Temperature derivative does not consider changes of electron temperature."); return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; } -void TwoTempPlasma::setContext(const Reaction& rxn, const Kinetics& kin) +void TwoTempPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) { // TwoTempPlasmaReaction is for a non-equilibrium plasma, and the reverse rate // cannot be calculated from the conventional thermochemistry. // @todo implement the reversible rate for non-equilibrium plasma if (rxn.reversible) { - throw InputFileError("TwoTempPlasma::setContext", rxn.input, + throw InputFileError("TwoTempPlasmaRate::setContext", rxn.input, "TwoTempPlasmaRate does not support reversible reactions"); } } diff --git a/test/kinetics/kineticsFromScratch3.cpp b/test/kinetics/kineticsFromScratch3.cpp index 7c912d564f1..06d489cc043 100644 --- a/test/kinetics/kineticsFromScratch3.cpp +++ b/test/kinetics/kineticsFromScratch3.cpp @@ -143,11 +143,11 @@ TEST_F(KineticsFromScratch3, add_plog_reaction) // - {P: 100.0 atm, A: 5.9632e+56, b: -11.529, Ea: 5.25996e+04} Composition reac = parseCompString("H2:1, O2:1"); Composition prod = parseCompString("OH:2"); - std::multimap rates { - { 0.01*101325, Arrhenius3(1.212400e+16, -0.5779, 10872.7 * 4184.0) }, - { 1.0*101325, Arrhenius3(4.910800e+31, -4.8507, 24772.8 * 4184.0) }, - { 10.0*101325, Arrhenius3(1.286600e+47, -9.0246, 39796.5 * 4184.0) }, - { 100.0*101325, Arrhenius3(5.963200e+56, -11.529, 52599.6 * 4184.0) } + std::multimap rates { + { 0.01*101325, ArrheniusRate(1.212400e+16, -0.5779, 10872.7 * 4184.0) }, + { 1.0*101325, ArrheniusRate(4.910800e+31, -4.8507, 24772.8 * 4184.0) }, + { 10.0*101325, ArrheniusRate(1.286600e+47, -9.0246, 39796.5 * 4184.0) }, + { 100.0*101325, ArrheniusRate(5.963200e+56, -11.529, 52599.6 * 4184.0) } }; auto R = make_shared(reac, prod, make_shared(rates)); @@ -159,11 +159,11 @@ TEST_F(KineticsFromScratch3, plog_invalid_rate) { Composition reac = parseCompString("H2:1, O2:1"); Composition prod = parseCompString("OH:2"); - std::multimap rates { - { 0.01*101325, Arrhenius3(1.2124e+16, -0.5779, 10872.7 * 4184.0) }, - { 10.0*101325, Arrhenius3(1e15, -1, 10000 * 4184.0) }, - { 10.0*101325, Arrhenius3(-2e20, -2.0, 20000 * 4184.0) }, - { 100.0*101325, Arrhenius3(5.9632e+56, -11.529, 52599.6 * 4184.0) } + std::multimap rates { + { 0.01*101325, ArrheniusRate(1.2124e+16, -0.5779, 10872.7 * 4184.0) }, + { 10.0*101325, ArrheniusRate(1e15, -1, 10000 * 4184.0) }, + { 10.0*101325, ArrheniusRate(-2e20, -2.0, 20000 * 4184.0) }, + { 100.0*101325, ArrheniusRate(5.9632e+56, -11.529, 52599.6 * 4184.0) } }; auto R = make_shared(reac, prod, make_shared(rates)); diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 5ab3ff1f8bb..9fd2d5356e0 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -659,7 +659,21 @@ TEST_F(ReactionToYaml, unconvertible2) EXPECT_THROW(params.applyUnits(), CanteraError); } -TEST_F(ReactionToYaml, BlowersMasel) +TEST_F(ReactionToYaml, unconvertible3) +{ + FalloffReaction3 R( + {{"H2", 1}, {"OH", 1}}, {{"H2O", 1}, {"H", 1}}, + TroeRate( + ArrheniusRate(1e5, -1.0, 12.5), ArrheniusRate(1e5, -1.0, 12.5), + {0.562, 91.0, 5836.0, 8552.0}), + ThirdBody()); + AnyMap params = R.parameters(); + UnitSystem U{"g", "cm", "mol"}; + params.setUnits(U); + EXPECT_THROW(params.applyUnits(), CanteraError); +} + +TEST_F(ReactionToYaml, BlowersMaselRate) { soln = newSolution("blowers-masel.yaml", "gas"); soln->thermo()->setState_TPY(1100, 0.1 * OneAtm, "O:0.01, H2:0.8, O2:0.19");