diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 3e381a2c95..eb74573aa2 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -149,6 +149,16 @@ 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 /*! @@ -203,213 +213,6 @@ class Arrhenius3 : public ArrheniusBase }; -//! Two temperature plasma reaction rate type depends on both -//! gas temperature and electron temperature. -/*! - * The form of the two temperature plasma reaction rate coefficient is similar to an - * Arrhenius reaction rate coefficient. The temperature exponent (b) is applied to - * the electron temperature instead. In addition, the exponential term with - * activation energy for electron is included. - * - * \f[ - * k_f = A T_e^b \exp (-E_a/RT) \exp (-E_{a,e}/RT_e) - * \f] - * - * Ref.: Kossyi, I. A., Kostinsky, A. Y., Matveyev, A. A., & Silakov, V. P. (1992). - * Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures. - * Plasma Sources Science and Technology, 1(3), 207. - * doi: 10.1088/0963-0252/1/3/011 - * - * @ingroup arrheniusGroup - */ -class TwoTempPlasma : public ArrheniusBase -{ -public: - TwoTempPlasma(); - - //! Constructor. - /*! - * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units - * depend on the reaction order and the dimensionality (surface or bulk). - * @param b Temperature exponent (non-dimensional) - * @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); - - virtual const std::string type() const override { - return "two-temperature-plasma"; - } - - virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; - - //! Evaluate reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double evalFromStruct(const TwoTempPlasmaData& shared_data) const { - // m_E4_R is the electron activation (in temperature units) - return m_A * std::exp(m_b * shared_data.logTe - - m_Ea_R * shared_data.recipT + - m_E4_R * (shared_data.electronTemp - shared_data.temperature) - * shared_data.recipTe * shared_data.recipT); - } - - //! Evaluate derivative of reaction rate with respect to temperature - //! divided by reaction rate - /*! - * This method does not consider changes of electron temperature. - * A corresponding warning is raised. - * @param shared_data data shared by all reactions of a given type - */ - double ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const; - - //! Return the electron activation energy *Ea* [J/kmol] - double activationElectronEnergy() const { - return m_E4_R * GasConstant; - } -}; - - -//! Blowers Masel reaction rate type depends on the enthalpy of reaction -/** - * The Blowers Masel approximation is written by Paul Blowers, - * Rich Masel (DOI: https://doi.org/10.1002/aic.690461015) to - * adjust the activation energy based on enthalpy change of a reaction: - * - * \f{eqnarray*}{ - * E_a &=& 0\; \text{if }\Delta H < -4E_0 \\ - * E_a &=& \Delta H\; \text{if }\Delta H > 4E_0 \\ - * E_a &=& \frac{(w + \Delta H / 2)(V_P - 2w + - * \Delta H)^2}{(V_P^2 - 4w^2 + (\Delta H)^2)}\; \text{Otherwise} - * \f} - * where - * \f[ - * V_P = \frac{2w (w + E_0)}{w - E_0}, - * \f] - * \f$ w \f$ is the average bond dissociation energy of the bond breaking - * and that being formed in the reaction. Since the expression is - * very insensitive to \f$ w \f$ for \f$ w >= 2 E_0 \f$, \f$ w \f$ - * can be approximated to an arbitrary high value like 1000 kJ/mol. - * - * After the activation energy is determined by Blowers-Masel approximation, - * it can be plugged into Arrhenius function to calculate the rate constant. - * \f[ - * k_f = A T^b \exp (-E_a/RT) - * \f] - * - * @ingroup arrheniusGroup - */ -class BlowersMasel : public ArrheniusBase -{ -public: - //! Default constructor. - BlowersMasel(); - - //! Constructor. - /*! - * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units - * depend on the reaction order and the dimensionality (surface or bulk). - * @param b Temperature exponent (non-dimensional) - * @param Ea0 Intrinsic activation energy in energy units [J/kmol] - * @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); - - virtual const std::string type() const override { - return "Blowers-Masel"; - } - - virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; - - //! Evaluate reaction rate - double evalRate(double logT, double recipT) const { - double Ea_R = effectiveActivationEnergy_R(m_deltaH_R); - return m_A * std::exp(m_b * logT - Ea_R * recipT); - } - - //! Update information specific to reaction - void updateFromStruct(const BlowersMaselData& shared_data) { - if (shared_data.ready) { - m_deltaH_R = 0.; - for (const auto& item : m_stoich_coeffs) { - m_deltaH_R += shared_data.partialMolarEnthalpies[item.first] * item.second; - } - m_deltaH_R /= GasConstant; - } - } - - //! Evaluate reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double evalFromStruct(const BlowersMaselData& shared_data) const { - double Ea_R = effectiveActivationEnergy_R(m_deltaH_R); - return m_A * std::exp(m_b * shared_data.logT - Ea_R * shared_data.recipT); - } - - //! divided by reaction rate - /*! - * This method does not consider potential changes due to a changed reaction - * enthalpy. A corresponding warning is raised. - * @param shared_data data shared by all reactions of a given type - */ - double ddTScaledFromStruct(const BlowersMaselData& shared_data) const; - -protected: - //! Return the effective activation energy (a function of the delta H of reaction) - //! divided by the gas constant (i.e. the activation temperature) [K] - //! @internal The enthalpy change of reaction is not an independent parameter - double effectiveActivationEnergy_R(double deltaH_R) const { - if (deltaH_R < -4 * m_Ea_R) { - return 0.; - } - if (deltaH_R > 4 * m_Ea_R) { - return deltaH_R; - } - // m_E4_R is the bond dissociation energy "w" (in temperature units) - double vp = 2 * m_E4_R * ((m_E4_R + m_Ea_R) / (m_E4_R - m_Ea_R)); // in Kelvin - double vp_2w_dH = (vp - 2 * m_E4_R + deltaH_R); // (Vp - 2 w + dH) - return (m_E4_R + deltaH_R / 2) * (vp_2w_dH * vp_2w_dH) / - (vp * vp - 4 * m_E4_R * m_E4_R + deltaH_R * deltaH_R); // in Kelvin - } - -public: - virtual double activationEnergy() const override { - return effectiveActivationEnergy_R(m_deltaH_R) * GasConstant; - } - - //! Return the bond dissociation energy *w* [J/kmol] - double bondEnergy() const { - return m_E4_R * GasConstant; - } - - //! Return current enthalpy change of reaction [J/kmol] - double deltaH() const { - return m_deltaH_R * GasConstant; - } - - //! Set current enthalpy change of reaction [J/kmol] - /*! - * @internal used for testing purposes only; note that this quantity is not an - * independent variable and will be overwritten during an update of the state. - * - * @warning This method is an experimental part of the %Cantera API and - * may be changed or removed without notice. - */ - void setDeltaH(double deltaH) { - m_deltaH_R = deltaH / GasConstant; - } - -protected: - //! Pairs of species indices and multiplers to calculate enthalpy change - std::vector> m_stoich_coeffs; - - double m_deltaH_R; //!< enthalpy change of reaction (in temperature units) -}; - - //! A class template for bulk phase reaction rate specifications template class BulkRate : public RateType @@ -456,8 +259,6 @@ class BulkRate : public RateType }; typedef BulkRate ArrheniusRate; -typedef BulkRate TwoTempPlasmaRate; -typedef BulkRate BlowersMaselRate; } diff --git a/include/cantera/kinetics/BlowersMaselRate.h b/include/cantera/kinetics/BlowersMaselRate.h new file mode 100644 index 0000000000..23cc818e27 --- /dev/null +++ b/include/cantera/kinetics/BlowersMaselRate.h @@ -0,0 +1,184 @@ +//! @file BlowersMaselRate.h Header for Blowers-Masel reaction rates + +// 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_BLOWERSMASELRATE_H +#define CT_BLOWERSMASELRATE_H + +#include "Arrhenius.h" + +namespace Cantera +{ + +//! Data container holding shared data specific to BlowersMaselRate +/** + * The data container `BlowersMaselData` holds precalculated data common to + * all `BlowersMaselRate` objects. + */ +struct BlowersMaselData : public ReactionData +{ + BlowersMaselData(); + + virtual void update(double T) override; + virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; + using ReactionData::update; + + virtual void resize(size_t nSpecies, size_t nReactions, size_t nPhases) override { + partialMolarEnthalpies.resize(nSpecies, 0.); + ready = true; + } + + bool ready; //!< boolean indicating whether vectors are accessible + double density; //!< used to determine if updates are needed + vector_fp partialMolarEnthalpies; //!< partial molar enthalpies + +protected: + int m_state_mf_number; //!< integer that is incremented when composition changes +}; + + +//! Blowers Masel reaction rate type depends on the enthalpy of reaction +/** + * The Blowers Masel approximation is written by Paul Blowers, + * Rich Masel (DOI: https://doi.org/10.1002/aic.690461015) to + * adjust the activation energy based on enthalpy change of a reaction: + * + * \f{eqnarray*}{ + * E_a &=& 0\; \text{if }\Delta H < -4E_0 \\ + * E_a &=& \Delta H\; \text{if }\Delta H > 4E_0 \\ + * E_a &=& \frac{(w + \Delta H / 2)(V_P - 2w + + * \Delta H)^2}{(V_P^2 - 4w^2 + (\Delta H)^2)}\; \text{Otherwise} + * \f} + * where + * \f[ + * V_P = \frac{2w (w + E_0)}{w - E_0}, + * \f] + * \f$ w \f$ is the average bond dissociation energy of the bond breaking + * and that being formed in the reaction. Since the expression is + * very insensitive to \f$ w \f$ for \f$ w >= 2 E_0 \f$, \f$ w \f$ + * can be approximated to an arbitrary high value like 1000 kJ/mol. + * + * After the activation energy is determined by Blowers-Masel approximation, + * it can be plugged into Arrhenius function to calculate the rate constant. + * \f[ + * k_f = A T^b \exp (-E_a/RT) + * \f] + * + * @ingroup arrheniusGroup + */ +class BlowersMasel : public ArrheniusBase +{ +public: + //! Default constructor. + BlowersMasel(); + + //! Constructor. + /*! + * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units + * depend on the reaction order and the dimensionality (surface or bulk). + * @param b Temperature exponent (non-dimensional) + * @param Ea0 Intrinsic activation energy in energy units [J/kmol] + * @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); + + virtual const std::string type() const override { + return "Blowers-Masel"; + } + + virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; + + //! Evaluate reaction rate + double evalRate(double logT, double recipT) const { + double Ea_R = effectiveActivationEnergy_R(m_deltaH_R); + return m_A * std::exp(m_b * logT - Ea_R * recipT); + } + + //! Update information specific to reaction + void updateFromStruct(const BlowersMaselData& shared_data) { + if (shared_data.ready) { + m_deltaH_R = 0.; + for (const auto& item : m_stoich_coeffs) { + m_deltaH_R += shared_data.partialMolarEnthalpies[item.first] * item.second; + } + m_deltaH_R /= GasConstant; + } + } + + //! Evaluate reaction rate + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const BlowersMaselData& shared_data) const { + double Ea_R = effectiveActivationEnergy_R(m_deltaH_R); + return m_A * std::exp(m_b * shared_data.logT - Ea_R * shared_data.recipT); + } + + //! Evaluate derivative of reaction rate with respect to temperature + //! divided by reaction rate + /*! + * This method does not consider potential changes due to a changed reaction + * enthalpy. A corresponding warning is raised. + * @param shared_data data shared by all reactions of a given type + */ + double ddTScaledFromStruct(const BlowersMaselData& shared_data) const; + +protected: + //! Return the effective activation energy (a function of the delta H of reaction) + //! divided by the gas constant (i.e. the activation temperature) [K] + //! @internal The enthalpy change of reaction is not an independent parameter + double effectiveActivationEnergy_R(double deltaH_R) const { + if (deltaH_R < -4 * m_Ea_R) { + return 0.; + } + if (deltaH_R > 4 * m_Ea_R) { + return deltaH_R; + } + // m_E4_R is the bond dissociation energy "w" (in temperature units) + double vp = 2 * m_E4_R * ((m_E4_R + m_Ea_R) / (m_E4_R - m_Ea_R)); // in Kelvin + double vp_2w_dH = (vp - 2 * m_E4_R + deltaH_R); // (Vp - 2 w + dH) + return (m_E4_R + deltaH_R / 2) * (vp_2w_dH * vp_2w_dH) / + (vp * vp - 4 * m_E4_R * m_E4_R + deltaH_R * deltaH_R); // in Kelvin + } + +public: + virtual double activationEnergy() const override { + return effectiveActivationEnergy_R(m_deltaH_R) * GasConstant; + } + + //! Return the bond dissociation energy *w* [J/kmol] + double bondEnergy() const { + return m_E4_R * GasConstant; + } + + //! Return current enthalpy change of reaction [J/kmol] + double deltaH() const { + return m_deltaH_R * GasConstant; + } + + //! Set current enthalpy change of reaction [J/kmol] + /*! + * @internal used for testing purposes only; note that this quantity is not an + * independent variable and will be overwritten during an update of the state. + * + * @warning This method is an experimental part of the %Cantera API and + * may be changed or removed without notice. + */ + void setDeltaH(double deltaH) { + m_deltaH_R = deltaH / GasConstant; + } + +protected: + //! Pairs of species indices and multiplers to calculate enthalpy change + std::vector> m_stoich_coeffs; + + 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 new file mode 100644 index 0000000000..37c66c4ad4 --- /dev/null +++ b/include/cantera/kinetics/ChebyshevRate.h @@ -0,0 +1,281 @@ +//! @file ChebyshevRate.h + +// 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_CHEBYSHEV_H +#define CT_CHEBYSHEV_H + +#include "cantera/kinetics/ReactionRate.h" +#include "cantera/kinetics/ReactionData.h" +#include "cantera/kinetics/MultiRate.h" +#include "cantera/base/Array.h" +#include "cantera/base/global.h" + +namespace Cantera +{ + +//! Data container holding shared data specific to ChebyshevRate +/** + * The data container `ChebyshevData` holds precalculated data common to + * all `ChebyshevRate` objects. + */ +struct ChebyshevData : public ReactionData +{ + ChebyshevData() : pressure(NAN), log10P(0.), m_pressure_buf(-1.) {} + + virtual void update(double T) override; + + virtual void update(double T, double P) override { + ReactionData::update(T); + pressure = P; + log10P = std::log10(P); + } + + virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; + + using ReactionData::update; + + //! Perturb pressure of data container + /** + * The method is used for the evaluation of numerical derivatives. + * @param deltaP relative pressure perturbation + */ + void perturbPressure(double deltaP); + + virtual void restore() override; + + virtual void invalidateCache() override { + ReactionData::invalidateCache(); + pressure = NAN; + } + + double pressure; //!< pressure + double log10P; //!< base 10 logarithm of pressure + +protected: + double m_pressure_buf; //!< buffered pressure +}; + +//! Pressure-dependent rate expression where the rate coefficient is expressed +//! as a bivariate Chebyshev polynomial in temperature and pressure. +/*! + * The rate constant can be written as: + * \f[ + * \log k(T,P) = \sum_{t=1}^{N_T} \sum_{p=1}^{N_P} \alpha_{tp} + * \phi_t(\tilde{T}) \phi_p(\tilde{P}) + * \f] + * where \f$\alpha_{tp}\f$ are the constants defining the rate, \f$\phi_n(x)\f$ + * is the Chebyshev polynomial of the first kind of degree *n* evaluated at + * *x*, and + * \f[ + * \tilde{T} \equiv \frac{2T^{-1} - T_\mathrm{min}^{-1} - T_\mathrm{max}^{-1}} + * {T_\mathrm{max}^{-1} - T_\mathrm{min}^{-1}} + * \f] + * \f[ + * \tilde{P} \equiv \frac{2 \log P - \log P_\mathrm{min} - \log P_\mathrm{max}} + * {\log P_\mathrm{max} - \log P_\mathrm{min}} + * \f] + * are reduced temperature and reduced pressures which map the ranges + * \f$ (T_\mathrm{min}, T_\mathrm{max}) \f$ and + * \f$ (P_\mathrm{min}, P_\mathrm{max}) \f$ to (-1, 1). + * + * A ChebyshevRate rate expression is specified in terms of the coefficient matrix + * \f$ \alpha \f$ and the temperature and pressure ranges. Note that the + * Chebyshev polynomials are not defined outside the interval (-1,1), and + * therefore extrapolation of rates outside the range of temperatures and + * pressures for which they are defined is strongly discouraged. + */ +class ChebyshevRate final : public ReactionRate +{ +public: + //! Default constructor. + ChebyshevRate() : m_log10P(NAN), m_rate_units(Units(0.)) {} + + //! Constructor directly from coefficient array + /*! + * @param Tmin Minimum temperature [K] + * @param Tmax Maximum temperature [K] + * @param Pmin Minimum pressure [Pa] + * @param Pmax Maximum pressure [Pa] + * @param coeffs Coefficient array dimensioned `nT` by `nP` where `nT` and + * `nP` are the number of temperatures and pressures used in the fit, + * respectively. + */ + ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs); + + ChebyshevRate(const AnyMap& node, const UnitStack& rate_units={}) + : ChebyshevRate() + { + setParameters(node, rate_units); + } + + unique_ptr newMultiRate() const { + return unique_ptr( + new MultiRate); + } + + const std::string type() const { return "Chebyshev"; } + + //! Perform object setup based on AnyMap node information + /*! + * @param node AnyMap containing rate information + * @param units Unit definitions specific to rate information + */ + void setParameters(const AnyMap& node, const UnitStack& units); + void getParameters(AnyMap& rateNode, const Units& rate_units) const { + // @todo: deprecate, as second argument is no longer needed + return getParameters(rateNode); + } + void getParameters(AnyMap& rateNode) const; + + //! Update information specific to reaction + /*! + * @param shared_data data shared by all reactions of a given type + */ + void updateFromStruct(const ChebyshevData& shared_data) { + if (shared_data.log10P != m_log10P) { + update_C(&shared_data.log10P); + } + } + + //! Evaluate reaction rate + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const ChebyshevData& shared_data) { + return updateRC(0., shared_data.recipT); + } + + //! Set up ChebyshevRate object + /*! + * @deprecated Deprecated in Cantera 2.6. Replaceable with + * @see setLimits() and @see setCoeffs(). + */ + void setup(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs); + + //! Set limits for ChebyshevRate object + /*! + * @param Tmin Minimum temperature [K] + * @param Tmax Maximum temperature [K] + * @param Pmin Minimum pressure [Pa] + * @param Pmax Maximum pressure [Pa] + */ + void setLimits(double Tmin, double Tmax, double Pmin, double Pmax); + + //! Update concentration-dependent parts of the rate coefficient. + //! @param c base-10 logarithm of the pressure in Pa + //! @deprecated To be removed after Cantera 2.6. Implementation will be moved to + //! the updateFromStruct() method. + void update_C(const double* c) { + m_log10P = c[0]; + double Pr = (2 * c[0] + PrNum_) * PrDen_; + double Cnm1 = Pr; + double Cn = 1; + double Cnp1; + for (size_t i = 0; i < m_coeffs.nRows(); i++) { + dotProd_[i] = m_coeffs(i, 0); + } + for (size_t j = 1; j < m_coeffs.nColumns(); j++) { + Cnp1 = 2 * Pr * Cn - Cnm1; + for (size_t i = 0; i < m_coeffs.nRows(); i++) { + dotProd_[i] += Cnp1 * m_coeffs(i, j); + } + Cnm1 = Cn; + Cn = Cnp1; + } + } + + /** + * Update the value the rate constant. + * + * This function returns the actual value of the rate constant. + * @deprecated To be removed after Cantera 2.6. Implementation will be moved to + * the evalFromStruct() method. + */ + double updateRC(double logT, double recipT) const { + double Tr = (2 * recipT + TrNum_) * TrDen_; + double Cnm1 = Tr; + double Cn = 1; + double Cnp1; + double logk = dotProd_[0]; + for (size_t i = 1; i < m_coeffs.nRows(); i++) { + Cnp1 = 2 * Tr * Cn - Cnm1; + logk += Cnp1 * dotProd_[i]; + Cnm1 = Cn; + Cn = Cnp1; + } + return std::pow(10, logk); + } + + //! Minimum valid temperature [K] + double Tmin() const { + return Tmin_; + } + + //! Maximum valid temperature [K] + double Tmax() const { + return Tmax_; + } + + //! Minimum valid pressure [Pa] + double Pmin() const { + return Pmin_; + } + + //! Maximum valid pressure [Pa] + double Pmax() const { + return Pmax_; + } + + //! Number of points in the pressure direction + size_t nPressure() const { + return m_coeffs.nColumns(); + } + + //! Number of points in the temperature direction + size_t nTemperature() const { + return m_coeffs.nRows(); + } + + //! Access the ChebyshevRate coefficients. + /*! + * \f$ \alpha_{t,p} = \mathrm{coeffs}[N_P*t + p] \f$ where + * \f$ 0 <= t < N_T \f$ and \f$ 0 <= p < N_P \f$. + * + * @deprecated To be removed after Cantera 2.6. Replaceable by @see data(). + */ + const vector_fp& coeffs() const { + warn_deprecated("ChebyshevRate::coeffs", "Deprecated in Cantera 2.6 " + "and to be removed thereafter; replaceable by data()."); + return chebCoeffs_; + } + + //! Access Chebyshev coefficients as 2-dimensional array with temperature and + //! pressure dimensions corresponding to rows and columns, respectively. + const Array2D& data() const { + return m_coeffs; + } + + //! Set the Chebyshev coefficients as 2-dimensional array. + void setData(const Array2D& coeffs); + +protected: + double m_log10P; //!< value detecting updates + double Tmin_, Tmax_; //!< valid temperature range + double Pmin_, Pmax_; //!< valid pressure range + double TrNum_, TrDen_; //!< terms appearing in the reduced temperature + double PrNum_, PrDen_; //!< terms appearing in the reduced pressure + + Array2D m_coeffs; //!<< coefficient array + vector_fp chebCoeffs_; //!< Chebyshev coefficients, length nP * nT + vector_fp dotProd_; //!< dot product of chebCoeffs with the reduced pressure polynomial + + Units m_rate_units; //!< Reaction rate units +}; + +} + +#endif diff --git a/include/cantera/kinetics/Custom.h b/include/cantera/kinetics/Custom.h index 7a7483e2be..6b5192ddfb 100644 --- a/include/cantera/kinetics/Custom.h +++ b/include/cantera/kinetics/Custom.h @@ -13,7 +13,7 @@ #include "cantera/base/ct_defs.h" #include "cantera/base/Units.h" -#include "cantera/kinetics/ReactionData.h" +#include "cantera/kinetics/Arrhenius.h" #include "ReactionRate.h" #include "MultiRate.h" diff --git a/include/cantera/kinetics/Falloff.h b/include/cantera/kinetics/Falloff.h index 799512cbe8..37efc67ae3 100644 --- a/include/cantera/kinetics/Falloff.h +++ b/include/cantera/kinetics/Falloff.h @@ -21,6 +21,55 @@ class AnyMap; * @ingroup chemkinetics */ + +//! Data container holding shared data specific to Falloff rates +/** + * The data container `FalloffData` holds precalculated data common to + * all Falloff related reaction rate classes. + */ +struct FalloffData : public ReactionData +{ + FalloffData(); + + virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; + + virtual void update(double T) override; + + virtual void update(double T, double M) override; + + using ReactionData::update; + + //! Perturb third-body concentration vector of data container + /** + * The method is used for the evaluation of numerical derivatives. + * @param deltaM relative third-body perturbation + */ + void perturbThirdBodies(double deltaM); + + virtual void restore() override; + + virtual void resize(size_t nSpecies, size_t nReactions, size_t nPhases) override { + conc_3b.resize(nReactions, NAN); + m_conc_3b_buf.resize(nReactions, NAN); + ready = true; + } + + virtual void invalidateCache() override { + ReactionData::invalidateCache(); + molar_density = NAN; + } + + bool ready; //!< boolean indicating whether vectors are accessible + double molar_density; //!< used to determine if updates are needed + vector_fp conc_3b; //!< vector of effective third-body concentrations + +protected: + int m_state_mf_number; //!< integer that is incremented when composition changes + bool m_perturbed; //!< boolean indicating whether 3-rd body values are perturbed + vector_fp m_conc_3b_buf; //!< buffered third-body concentrations +}; + + /** * Base class for falloff rate calculators. Each instance of a subclass of FalloffRate * calculates the falloff reaction rate based on specific implementations of the diff --git a/include/cantera/kinetics/InterfaceRate.h b/include/cantera/kinetics/InterfaceRate.h index 43a93387cd..7c21b328ed 100644 --- a/include/cantera/kinetics/InterfaceRate.h +++ b/include/cantera/kinetics/InterfaceRate.h @@ -10,7 +10,7 @@ #define CT_INTERFACERATE_H #include "cantera/base/global.h" -#include "cantera/kinetics/Arrhenius.h" +#include "cantera/kinetics/BlowersMaselRate.h" #include "MultiRate.h" namespace Cantera @@ -27,6 +27,49 @@ class AnyMap; * @ingroup chemkinetics */ + +//! Data container holding shared data for reaction rate specification with interfaces +/** + * The data container InterfaceData holds precalculated data common to + * InterfaceRate and StickingRate objects. + * + * The data container inherits from BlowersMaselData, where density is used to + * hold the site density [kmol/m^2]. + */ +struct InterfaceData : public BlowersMaselData +{ + InterfaceData(); + + virtual bool update(const ThermoPhase& bulk, const Kinetics& kin) override; + + virtual void update(double T) override; + + virtual void update(double T, const vector_fp& values) override; + + using BlowersMaselData::update; + + virtual void perturbTemperature(double deltaT); + + virtual void resize(size_t nSpecies, size_t nReactions, size_t nPhases) override { + coverages.resize(nSpecies, 0.); + logCoverages.resize(nSpecies, 0.); + partialMolarEnthalpies.resize(nSpecies, 0.); + electricPotentials.resize(nPhases, 0.); + standardChemPotentials.resize(nSpecies, 0.); + standardConcentrations.resize(nSpecies, 0.); + ready = true; + } + + double sqrtT; //!< square root of temperature + + vector_fp coverages; //!< surface coverages + vector_fp logCoverages; //!< logarithm of surface coverages + vector_fp electricPotentials; //!< electric potentials of phases + vector_fp standardChemPotentials; //!< standard state chemical potentials + vector_fp standardConcentrations; //!< standard state concentrations +}; + + //! Base class for rate parameterizations that involve interfaces /** * Rate expressions defined for interfaces may include coverage dependent terms, diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 74b8a3cf5b..44da8ef3e5 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -13,7 +13,6 @@ #include "StoichManager.h" #include "cantera/base/ValueCache.h" -#include "cantera/kinetics/ReactionFactory.h" namespace Cantera { diff --git a/include/cantera/kinetics/PlogRate.h b/include/cantera/kinetics/PlogRate.h new file mode 100644 index 0000000000..128a979f39 --- /dev/null +++ b/include/cantera/kinetics/PlogRate.h @@ -0,0 +1,251 @@ +//! @file PlogRate.h + +// 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_PLOGRATE_H +#define CT_PLOGRATE_H + +#include "cantera/kinetics/Arrhenius.h" + +namespace Cantera +{ + +class Arrhenius2; + +//! Data container holding shared data specific to PlogRate +/** + * The data container `PlogData` holds precalculated data common to + * all `PlogRate` objects. + */ +struct PlogData : public ReactionData +{ + PlogData() : pressure(NAN), logP(0.), m_pressure_buf(-1.) {} + + virtual void update(double T) override; + + virtual void update(double T, double P) override { + ReactionData::update(T); + pressure = P; + logP = std::log(P); + } + + virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; + + using ReactionData::update; + + //! Perturb pressure of data container + /** + * The method is used for the evaluation of numerical derivatives. + * @param deltaP relative pressure perturbation + */ + void perturbPressure(double deltaP); + + virtual void restore() override; + + virtual void invalidateCache() override { + ReactionData::invalidateCache(); + pressure = NAN; + } + + double pressure; //!< pressure + double logP; //!< logarithm of pressure + +protected: + double m_pressure_buf; //!< buffered pressure +}; + + +//! Pressure-dependent reaction rate expressed by logarithmically interpolating +//! between Arrhenius rate expressions at various pressures. +/*! + * Given two rate expressions at two specific pressures: + * + * * \f$ P_1: k_1(T) = A_1 T^{b_1} e^{-E_1 / RT} \f$ + * * \f$ P_2: k_2(T) = A_2 T^{b_2} e^{-E_2 / RT} \f$ + * + * The rate at an intermediate pressure \f$ P_1 < P < P_2 \f$ is computed as + * \f[ + * \log k(T,P) = \log k_1(T) + \bigl(\log k_2(T) - \log k_1(T)\bigr) + * \frac{\log P - \log P_1}{\log P_2 - \log P_1} + * \f] + * Multiple rate expressions may be given at the same pressure, in which case + * the rate used in the interpolation formula is the sum of all the rates given + * at that pressure. For pressures outside the given range, the rate expression + * at the nearest pressure is used. + */ +class PlogRate final : public ReactionRate +{ +public: + //! Default constructor. + PlogRate(); + + //! Constructor from Arrhenius rate expressions at a set of pressures + explicit PlogRate(const std::multimap& rates); + + //! Constructor using legacy Arrhenius2 framework + explicit PlogRate(const std::multimap& rates); + + PlogRate(const AnyMap& node, const UnitStack& rate_units={}) : PlogRate() { + setParameters(node, rate_units); + } + + unique_ptr newMultiRate() const { + return unique_ptr(new MultiRate); + } + + //! Identifier of reaction rate type + const std::string type() const { return "pressure-dependent-Arrhenius"; } + + //! Perform object setup based on AnyMap node information + /*! + * @param node AnyMap containing rate information + * @param units Unit definitions specific to rate information + */ + 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)); + } + + //! Update information specific to reaction + /*! + * @param shared_data data shared by all reactions of a given type + */ + void updateFromStruct(const PlogData& shared_data) { + if (shared_data.logP != logP_) { + update_C(&shared_data.logP); + } + } + + //! Evaluate reaction rate + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const PlogData& shared_data) { + return updateRC(shared_data.logT, shared_data.recipT); + } + + //! Set up Plog object + /*! + * @deprecated Deprecated in Cantera 2.6. Replaced by setRates. + */ + void setup(const std::multimap& rates); + + //! Set up Plog object + void setRates(const std::multimap& rates); + + //! Update concentration-dependent parts of the rate coefficient. + //! @param c natural log of the pressure in Pa + //! @deprecated To be removed after Cantera 2.6. Implementation will be moved to + //! the updateFromStruct() method. + void update_C(const double* c) { + logP_ = c[0]; + if (logP_ > logP1_ && logP_ < logP2_) { + return; + } + + auto iter = pressures_.upper_bound(c[0]); + AssertThrowMsg(iter != pressures_.end(), "PlogRate::update_C", + "Pressure out of range: {}", logP_); + AssertThrowMsg(iter != pressures_.begin(), "PlogRate::update_C", + "Pressure out of range: {}", logP_); + + // upper interpolation pressure + logP2_ = iter->first; + ihigh1_ = iter->second.first; + ihigh2_ = iter->second.second; + + // lower interpolation pressure + logP1_ = (--iter)->first; + ilow1_ = iter->second.first; + ilow2_ = iter->second.second; + + rDeltaP_ = 1.0 / (logP2_ - logP1_); + } + + /** + * Update the value the rate constant. + * + * This function returns the actual value of the rate constant. + * @deprecated To be removed after Cantera 2.6. Implementation will be moved to + * the evalFromStruct() method. + */ + double updateRC(double logT, double recipT) const { + double log_k1, log_k2; + if (ilow1_ == ilow2_) { + log_k1 = rates_[ilow1_].evalLog(logT, recipT); + } else { + double k = 1e-300; // non-zero to make log(k) finite + for (size_t i = ilow1_; i < ilow2_; i++) { + k += rates_[i].evalRate(logT, recipT); + } + log_k1 = std::log(k); + } + + if (ihigh1_ == ihigh2_) { + log_k2 = rates_[ihigh1_].evalLog(logT, recipT); + } else { + double k = 1e-300; // non-zero to make log(k) finite + for (size_t i = ihigh1_; i < ihigh2_; i++) { + k += rates_[i].evalRate(logT, recipT); + } + log_k2 = std::log(k); + } + + return std::exp(log_k1 + (log_k2-log_k1) * (logP_-logP1_) * rDeltaP_); + } + + //! Check to make sure that the rate expression is finite over a range of + //! temperatures at each interpolation pressure. This is potentially an + //! issue when one of the Arrhenius expressions at a particular pressure + //! has a negative pre-exponential factor. + void validate(const std::string& equation, const Kinetics& kin) { + validate(equation); + } + + void validate(const std::string& equation); + + //! Return the pressures and Arrhenius expressions which comprise this + //! reaction. + /*! + * @deprecated Behavior to change after Cantera 2.6. + * @see getRates for new behavior. + */ + std::vector> rates() const; + + //! Return the pressures and Arrhenius expressions which comprise this + //! reaction. + 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_; + + double logP_; //!< log(p) at the current state + double logP1_, logP2_; //!< log(p) at the lower / upper pressure reference + + //! Indices to the ranges within rates_ for the lower / upper pressure, such + //! that rates_[ilow1_] through rates_[ilow2_] (inclusive) are the rates + //! expressions which are combined to form the rate at the lower reference + //! pressure. + size_t ilow1_, ilow2_, ihigh1_, ihigh2_; + + double rDeltaP_; //!< reciprocal of (logP2 - logP1) +}; + +} +#endif \ No newline at end of file diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 59a94ccd5b..4a644e9394 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -12,6 +12,7 @@ #include "cantera/kinetics/ReactionRate.h" #include "cantera/kinetics/RxnRates.h" #include "cantera/base/Units.h" +#include "ChebyshevRate.h" #include "InterfaceRate.h" #include "Custom.h" @@ -554,6 +555,26 @@ typedef InterfaceReaction2 InterfaceReaction; typedef ElectrochemicalReaction2 ElectrochemicalReaction; #endif +//! Create a new empty Reaction object +/*! + * @param type string identifying type of reaction. + * @deprecated To be removed after Cantera 2.6. Only used for legacy reaction types. + */ +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); //! Create Reaction objects for all `` nodes in an XML document. //! diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h index ecc4a1c0d6..9cbe1ae861 100644 --- a/include/cantera/kinetics/ReactionData.h +++ b/include/cantera/kinetics/ReactionData.h @@ -9,6 +9,7 @@ #define CT_REACTIONDATA_H #include "cantera/base/ct_defs.h" +#include "cantera/base/ctexceptions.h" namespace Cantera { @@ -43,7 +44,10 @@ struct ReactionData * This method allows for testing of a reaction rate expression outside of * Kinetics reaction rate evaluators. */ - virtual void update(double T, double extra); + virtual void update(double T, double extra) { + throw NotImplementedError("ReactionData::update", + "ReactionData type does not use extra scalar argument."); + } //! Update data container based on temperature *T* and a vector parameter *extra* /** @@ -54,7 +58,10 @@ struct ReactionData * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - virtual void update(double T, const vector_fp& extra); + virtual void update(double T, const vector_fp& extra) { + throw NotImplementedError("ReactionData::update", + "ReactionData type does not use extra vector argument."); + } //! Update data container based on thermodynamic phase state /** @@ -69,10 +76,24 @@ struct ReactionData * The method is used for the evaluation of numerical derivatives. * @param deltaT relative temperature perturbation */ - void perturbTemperature(double deltaT); + void perturbTemperature(double deltaT) { + if (m_temperature_buf > 0.) { + throw CanteraError("ReactionData::perturbTemperature", + "Cannot apply another perturbation as state is already perturbed."); + } + m_temperature_buf = temperature; + ReactionData::update(temperature * (1. + deltaT)); + } //! Restore data container after a perturbation - virtual void restore(); + virtual void restore() { + // only restore if there is a valid buffered value + if (m_temperature_buf < 0.) { + return; + } + ReactionData::update(m_temperature_buf); + m_temperature_buf = -1.; + } //! Update number of species, reactions and phases virtual void resize(size_t nSpecies, size_t nReactions, size_t nPhases) {} @@ -93,253 +114,6 @@ struct ReactionData double m_temperature_buf; //!< buffered temperature }; - -//! 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; -}; - - -//! Data container holding shared data specific to TwoTempPlasmaRate -/** - * The data container `TwoTempPlasmaData` holds precalculated data common to - * all `TwoTempPlasmaRate` objects. - */ -struct TwoTempPlasmaData : public ReactionData -{ - TwoTempPlasmaData() : electronTemp(1.), logTe(0.), recipTe(1.) {} - - virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; - - virtual void update(double T) override; - - virtual void update(double T, double Te) override; - - using ReactionData::update; - - virtual void updateTe(double Te); - - virtual void invalidateCache() override { - ReactionData::invalidateCache(); - electronTemp = NAN; - } - - double electronTemp; //!< electron temperature - double logTe; //!< logarithm of electron temperature - double recipTe; //!< inverse of electron temperature -}; - - -//! Data container holding shared data specific to BlowersMaselRate -/** - * The data container `BlowersMaselData` holds precalculated data common to - * all `BlowersMaselRate` objects. - */ -struct BlowersMaselData : public ReactionData -{ - BlowersMaselData(); - - virtual void update(double T) override; - - virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; - - using ReactionData::update; - - virtual void resize(size_t nSpecies, size_t nReactions, size_t nPhases) override { - partialMolarEnthalpies.resize(nSpecies, 0.); - ready = true; - } - - bool ready; //!< boolean indicating whether vectors are accessible - double density; //!< used to determine if updates are needed - vector_fp partialMolarEnthalpies; //!< partial molar enthalpies - -protected: - int m_state_mf_number; //!< integer that is incremented when composition changes -}; - - -//! Data container holding shared data specific to Falloff rates -/** - * The data container `FalloffData` holds precalculated data common to - * all Falloff related reaction rate classes. - */ -struct FalloffData : public ReactionData -{ - FalloffData(); - - virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; - - virtual void update(double T) override; - - virtual void update(double T, double M) override; - - using ReactionData::update; - - //! Perturb third-body concentration vector of data container - /** - * The method is used for the evaluation of numerical derivatives. - * @param deltaM relative third-body perturbation - */ - void perturbThirdBodies(double deltaM); - - virtual void restore() override; - - virtual void resize(size_t nSpecies, size_t nReactions, size_t nPhases) override { - conc_3b.resize(nReactions, NAN); - m_conc_3b_buf.resize(nReactions, NAN); - ready = true; - } - - virtual void invalidateCache() override { - ReactionData::invalidateCache(); - molar_density = NAN; - } - - bool ready; //!< boolean indicating whether vectors are accessible - double molar_density; //!< used to determine if updates are needed - vector_fp conc_3b; //!< vector of effective third-body concentrations - -protected: - int m_state_mf_number; //!< integer that is incremented when composition changes - bool m_perturbed; //!< boolean indicating whether 3-rd body values are perturbed - vector_fp m_conc_3b_buf; //!< buffered third-body concentrations -}; - - -//! Data container holding shared data specific to PlogRate -/** - * The data container `PlogData` holds precalculated data common to - * all `PlogRate` objects. - */ -struct PlogData : public ReactionData -{ - PlogData() : pressure(NAN), logP(0.), m_pressure_buf(-1.) {} - - virtual void update(double T) override; - - virtual void update(double T, double P) override { - ReactionData::update(T); - pressure = P; - logP = std::log(P); - } - - virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; - - using ReactionData::update; - - //! Perturb pressure of data container - /** - * The method is used for the evaluation of numerical derivatives. - * @param deltaP relative pressure perturbation - */ - void perturbPressure(double deltaP); - - virtual void restore() override; - - virtual void invalidateCache() override { - ReactionData::invalidateCache(); - pressure = NAN; - } - - double pressure; //!< pressure - double logP; //!< logarithm of pressure - -protected: - double m_pressure_buf; //!< buffered pressure -}; - - -//! Data container holding shared data specific to ChebyshevRate -/** - * The data container `ChebyshevData` holds precalculated data common to - * all `ChebyshevRate` objects. - */ -struct ChebyshevData : public ReactionData -{ - ChebyshevData() : pressure(NAN), log10P(0.), m_pressure_buf(-1.) {} - - virtual void update(double T) override; - - virtual void update(double T, double P) override { - ReactionData::update(T); - pressure = P; - log10P = std::log10(P); - } - - virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; - - using ReactionData::update; - - //! Perturb pressure of data container - /** - * The method is used for the evaluation of numerical derivatives. - * @param deltaP relative pressure perturbation - */ - void perturbPressure(double deltaP); - - virtual void restore() override; - - virtual void invalidateCache() override { - ReactionData::invalidateCache(); - pressure = NAN; - } - - double pressure; //!< pressure - double log10P; //!< base 10 logarithm of pressure - -protected: - double m_pressure_buf; //!< buffered pressure -}; - - -//! Data container holding shared data for reaction rate specification with interfaces -/** - * The data container InterfaceData holds precalculated data common to - * InterfaceRate and StickingRate objects. - * - * The data container inherits from BlowersMaselData, where density is used to - * hold the site density [kmol/m^2]. - */ -struct InterfaceData : public BlowersMaselData -{ - InterfaceData(); - - virtual bool update(const ThermoPhase& bulk, const Kinetics& kin) override; - - virtual void update(double T) override; - - virtual void update(double T, const vector_fp& values) override; - - using BlowersMaselData::update; - - virtual void perturbTemperature(double deltaT); - - virtual void resize(size_t nSpecies, size_t nReactions, size_t nPhases) override { - coverages.resize(nSpecies, 0.); - logCoverages.resize(nSpecies, 0.); - partialMolarEnthalpies.resize(nSpecies, 0.); - electricPotentials.resize(nPhases, 0.); - standardChemPotentials.resize(nSpecies, 0.); - standardConcentrations.resize(nSpecies, 0.); - ready = true; - } - - double sqrtT; //!< square root of temperature - - vector_fp coverages; //!< surface coverages - vector_fp logCoverages; //!< logarithm of surface coverages - vector_fp electricPotentials; //!< electric potentials of phases - vector_fp standardChemPotentials; //!< standard state chemical potentials - vector_fp standardConcentrations; //!< standard state concentrations -}; - } #endif diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 80c4fb0a07..816ca97bb7 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -1,8 +1,7 @@ /** * @file RxnRates.h - * - * @todo at least part of the content of this file should be transferred - * to ReactionRate.h once the old XML interface is removed after Cantera 2.6 + * @deprecated To be removed after Cantera 2.6. See class Cantera::ReactionRate and + * derived classes for new reaction rate handlers. */ // This file is part of Cantera. See License.txt in the top-level directory or @@ -13,6 +12,7 @@ #include "cantera/kinetics/reaction_defs.h" #include "cantera/kinetics/Arrhenius.h" +#include "cantera/kinetics/PlogRate.h" #include "MultiRate.h" #include "cantera/base/Array.h" #include "cantera/base/ctexceptions.h" @@ -35,6 +35,8 @@ class Func1; * \f[ * k_f = A T^b \exp (-E/RT) * \f] + * + * @deprecated To be removed after Cantera 2.6. See Arrhenius3 / ArrheniusRate. */ class Arrhenius2 final : public Arrhenius3 { @@ -126,6 +128,8 @@ class Arrhenius2 final : public Arrhenius3 * \f] * 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$. + * + * @deprecated To be removed after Cantera 2.6. See InterfaceRate and StickingRate. */ class SurfaceArrhenius { @@ -203,411 +207,8 @@ typedef Arrhenius3 Arrhenius; typedef Arrhenius2 Arrhenius; #endif - -//! Pressure-dependent reaction rate expressed by logarithmically interpolating -//! between Arrhenius rate expressions at various pressures. -/*! - * Given two rate expressions at two specific pressures: - * - * * \f$ P_1: k_1(T) = A_1 T^{b_1} e^{-E_1 / RT} \f$ - * * \f$ P_2: k_2(T) = A_2 T^{b_2} e^{-E_2 / RT} \f$ - * - * The rate at an intermediate pressure \f$ P_1 < P < P_2 \f$ is computed as - * \f[ - * \log k(T,P) = \log k_1(T) + \bigl(\log k_2(T) - \log k_1(T)\bigr) - * \frac{\log P - \log P_1}{\log P_2 - \log P_1} - * \f] - * Multiple rate expressions may be given at the same pressure, in which case - * the rate used in the interpolation formula is the sum of all the rates given - * at that pressure. For pressures outside the given range, the rate expression - * at the nearest pressure is used. - */ -class PlogRate final : public ReactionRate -{ -public: - //! Default constructor. - PlogRate(); - - //! Constructor from Arrhenius rate expressions at a set of pressures - explicit PlogRate(const std::multimap& rates); - - //! Constructor using legacy Arrhenius2 framework - explicit PlogRate(const std::multimap& rates); - - PlogRate(const AnyMap& node, const UnitStack& rate_units={}) : PlogRate() { - setParameters(node, rate_units); - } - - unique_ptr newMultiRate() const { - return unique_ptr(new MultiRate); - } - - //! Identifier of reaction rate type - const std::string type() const { return "pressure-dependent-Arrhenius"; } - - //! Perform object setup based on AnyMap node information - /*! - * @param node AnyMap containing rate information - * @param units Unit definitions specific to rate information - */ - 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)); - } - - //! Update information specific to reaction - /*! - * @param shared_data data shared by all reactions of a given type - */ - void updateFromStruct(const PlogData& shared_data) { - if (shared_data.logP != logP_) { - update_C(&shared_data.logP); - } - } - - //! Evaluate reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double evalFromStruct(const PlogData& shared_data) { - return updateRC(shared_data.logT, shared_data.recipT); - } - - //! Set up Plog object - /*! - * @deprecated Deprecated in Cantera 2.6. Replaced by setRates. - */ - void setup(const std::multimap& rates); - - //! Set up Plog object - void setRates(const std::multimap& rates); - - //! Update concentration-dependent parts of the rate coefficient. - //! @param c natural log of the pressure in Pa - void update_C(const doublereal* c) { - logP_ = c[0]; - if (logP_ > logP1_ && logP_ < logP2_) { - return; - } - - auto iter = pressures_.upper_bound(c[0]); - AssertThrowMsg(iter != pressures_.end(), "PlogRate::update_C", - "Pressure out of range: {}", logP_); - AssertThrowMsg(iter != pressures_.begin(), "PlogRate::update_C", - "Pressure out of range: {}", logP_); - - // upper interpolation pressure - logP2_ = iter->first; - ihigh1_ = iter->second.first; - ihigh2_ = iter->second.second; - - // lower interpolation pressure - logP1_ = (--iter)->first; - ilow1_ = iter->second.first; - ilow2_ = iter->second.second; - - rDeltaP_ = 1.0 / (logP2_ - logP1_); - } - - /** - * Update the value the rate constant. - * - * This function returns the actual value of the rate constant. - */ - doublereal updateRC(doublereal logT, doublereal recipT) const { - double log_k1, log_k2; - if (ilow1_ == ilow2_) { - log_k1 = rates_[ilow1_].evalLog(logT, recipT); - } else { - double k = 1e-300; // non-zero to make log(k) finite - for (size_t i = ilow1_; i < ilow2_; i++) { - k += rates_[i].evalRate(logT, recipT); - } - log_k1 = std::log(k); - } - - if (ihigh1_ == ihigh2_) { - log_k2 = rates_[ihigh1_].evalLog(logT, recipT); - } else { - double k = 1e-300; // non-zero to make log(k) finite - for (size_t i = ihigh1_; i < ihigh2_; i++) { - k += rates_[i].evalRate(logT, recipT); - } - log_k2 = std::log(k); - } - - return std::exp(log_k1 + (log_k2-log_k1) * (logP_-logP1_) * rDeltaP_); - } - - //! Check to make sure that the rate expression is finite over a range of - //! temperatures at each interpolation pressure. This is potentially an - //! issue when one of the Arrhenius expressions at a particular pressure - //! has a negative pre-exponential factor. - void validate(const std::string& equation, const Kinetics& kin) { - validate(equation); - } - - void validate(const std::string& equation); - - //! Return the pressures and Arrhenius expressions which comprise this - //! reaction. - /*! - * @deprecated Behavior to change after Cantera 2.6. - * @see getRates for new behavior. - */ - std::vector> rates() const; - - //! Return the pressures and Arrhenius expressions which comprise this - //! reaction. - 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_; - - double logP_; //!< log(p) at the current state - double logP1_, logP2_; //!< log(p) at the lower / upper pressure reference - - //! Indices to the ranges within rates_ for the lower / upper pressure, such - //! that rates_[ilow1_] through rates_[ilow2_] (inclusive) are the rates - //! expressions which are combined to form the rate at the lower reference - //! pressure. - size_t ilow1_, ilow2_, ihigh1_, ihigh2_; - - double rDeltaP_; //!< reciprocal of (logP2 - logP1) -}; - typedef PlogRate Plog; -//! Pressure-dependent rate expression where the rate coefficient is expressed -//! as a bivariate Chebyshev polynomial in temperature and pressure. -/*! - * The rate constant can be written as: - * \f[ - * \log k(T,P) = \sum_{t=1}^{N_T} \sum_{p=1}^{N_P} \alpha_{tp} - * \phi_t(\tilde{T}) \phi_p(\tilde{P}) - * \f] - * where \f$\alpha_{tp}\f$ are the constants defining the rate, \f$\phi_n(x)\f$ - * is the Chebyshev polynomial of the first kind of degree *n* evaluated at - * *x*, and - * \f[ - * \tilde{T} \equiv \frac{2T^{-1} - T_\mathrm{min}^{-1} - T_\mathrm{max}^{-1}} - * {T_\mathrm{max}^{-1} - T_\mathrm{min}^{-1}} - * \f] - * \f[ - * \tilde{P} \equiv \frac{2 \log P - \log P_\mathrm{min} - \log P_\mathrm{max}} - * {\log P_\mathrm{max} - \log P_\mathrm{min}} - * \f] - * are reduced temperature and reduced pressures which map the ranges - * \f$ (T_\mathrm{min}, T_\mathrm{max}) \f$ and - * \f$ (P_\mathrm{min}, P_\mathrm{max}) \f$ to (-1, 1). - * - * A ChebyshevRate rate expression is specified in terms of the coefficient matrix - * \f$ \alpha \f$ and the temperature and pressure ranges. Note that the - * Chebyshev polynomials are not defined outside the interval (-1,1), and - * therefore extrapolation of rates outside the range of temperatures and - * pressures for which they are defined is strongly discouraged. - */ -class ChebyshevRate final : public ReactionRate -{ -public: - //! Default constructor. - ChebyshevRate() : m_log10P(NAN), m_rate_units(Units(0.)) {} - - //! Constructor directly from coefficient array - /*! - * @param Tmin Minimum temperature [K] - * @param Tmax Maximum temperature [K] - * @param Pmin Minimum pressure [Pa] - * @param Pmax Maximum pressure [Pa] - * @param coeffs Coefficient array dimensioned `nT` by `nP` where `nT` and - * `nP` are the number of temperatures and pressures used in the fit, - * respectively. - */ - ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, - const Array2D& coeffs); - - ChebyshevRate(const AnyMap& node, const UnitStack& rate_units={}) - : ChebyshevRate() - { - setParameters(node, rate_units); - } - - unique_ptr newMultiRate() const { - return unique_ptr( - new MultiRate); - } - - const std::string type() const { return "Chebyshev"; } - - //! Perform object setup based on AnyMap node information - /*! - * @param node AnyMap containing rate information - * @param units Unit definitions specific to rate information - */ - void setParameters(const AnyMap& node, const UnitStack& units); - void getParameters(AnyMap& rateNode, const Units& rate_units) const { - // @todo: deprecate, as second argument is no longer needed - return getParameters(rateNode); - } - void getParameters(AnyMap& rateNode) const; - - //! Update information specific to reaction - /*! - * @param shared_data data shared by all reactions of a given type - */ - void updateFromStruct(const ChebyshevData& shared_data) { - if (shared_data.log10P != m_log10P) { - update_C(&shared_data.log10P); - } - } - - //! Evaluate reaction rate - /*! - * @param shared_data data shared by all reactions of a given type - */ - double evalFromStruct(const ChebyshevData& shared_data) { - return updateRC(0., shared_data.recipT); - } - - //! Set up ChebyshevRate object - /*! - * @deprecated Deprecated in Cantera 2.6. Replaceable with - * @see setLimits() and @see setCoeffs(). - */ - void setup(double Tmin, double Tmax, double Pmin, double Pmax, - const Array2D& coeffs); - - //! Set limits for ChebyshevRate object - /*! - * @param Tmin Minimum temperature [K] - * @param Tmax Maximum temperature [K] - * @param Pmin Minimum pressure [Pa] - * @param Pmax Maximum pressure [Pa] - */ - void setLimits(double Tmin, double Tmax, double Pmin, double Pmax); - - //! Update concentration-dependent parts of the rate coefficient. - //! @param c base-10 logarithm of the pressure in Pa - void update_C(const double* c) { - m_log10P = c[0]; - double Pr = (2 * c[0] + PrNum_) * PrDen_; - double Cnm1 = Pr; - double Cn = 1; - double Cnp1; - for (size_t i = 0; i < m_coeffs.nRows(); i++) { - dotProd_[i] = m_coeffs(i, 0); - } - for (size_t j = 1; j < m_coeffs.nColumns(); j++) { - Cnp1 = 2 * Pr * Cn - Cnm1; - for (size_t i = 0; i < m_coeffs.nRows(); i++) { - dotProd_[i] += Cnp1 * m_coeffs(i, j); - } - Cnm1 = Cn; - Cn = Cnp1; - } - } - - /** - * Update the value the rate constant. - * - * This function returns the actual value of the rate constant. - */ - double updateRC(double logT, double recipT) const { - double Tr = (2 * recipT + TrNum_) * TrDen_; - double Cnm1 = Tr; - double Cn = 1; - double Cnp1; - double logk = dotProd_[0]; - for (size_t i = 1; i < m_coeffs.nRows(); i++) { - Cnp1 = 2 * Tr * Cn - Cnm1; - logk += Cnp1 * dotProd_[i]; - Cnm1 = Cn; - Cn = Cnp1; - } - return std::pow(10, logk); - } - - //! Minimum valid temperature [K] - double Tmin() const { - return Tmin_; - } - - //! Maximum valid temperature [K] - double Tmax() const { - return Tmax_; - } - - //! Minimum valid pressure [Pa] - double Pmin() const { - return Pmin_; - } - - //! Maximum valid pressure [Pa] - double Pmax() const { - return Pmax_; - } - - //! Number of points in the pressure direction - size_t nPressure() const { - return m_coeffs.nColumns(); - } - - //! Number of points in the temperature direction - size_t nTemperature() const { - return m_coeffs.nRows(); - } - - //! Access the ChebyshevRate coefficients. - /*! - * \f$ \alpha_{t,p} = \mathrm{coeffs}[N_P*t + p] \f$ where - * \f$ 0 <= t < N_T \f$ and \f$ 0 <= p < N_P \f$. - * - * @deprecated To be removed after Cantera 2.6. Replaceable by @see data(). - */ - const vector_fp& coeffs() const { - warn_deprecated("ChebyshevRate::coeffs", "Deprecated in Cantera 2.6 " - "and to be removed thereafter; replaceable by data()."); - return chebCoeffs_; - } - - //! Access Chebyshev coefficients as 2-dimensional array with temperature and - //! pressure dimensions corresponding to rows and columns, respectively. - const Array2D& data() const { - return m_coeffs; - } - - //! Set the Chebyshev coefficients as 2-dimensional array. - void setData(const Array2D& coeffs); - -protected: - double m_log10P; //!< value detecting updates - double Tmin_, Tmax_; //!< valid temperature range - double Pmin_, Pmax_; //!< valid pressure range - double TrNum_, TrDen_; //!< terms appearing in the reduced temperature - double PrNum_, PrDen_; //!< terms appearing in the reduced pressure - - Array2D m_coeffs; //!<< coefficient array - vector_fp chebCoeffs_; //!< Chebyshev coefficients, length nP * nT - vector_fp dotProd_; //!< dot product of chebCoeffs with the reduced pressure polynomial - - Units m_rate_units; //!< Reaction rate units -}; - } #endif diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h new file mode 100644 index 0000000000..57cb567e41 --- /dev/null +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -0,0 +1,113 @@ +//! @file TwoTempPlasmaRate.h Header for plasma reaction rates parameterized by two +//! temperatures (gas and electron). + +// 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_TWOTEMPPLASMARATE_H +#define CT_TWOTEMPPLASMARATE_H + +#include "Arrhenius.h" + +namespace Cantera +{ + +//! Data container holding shared data specific to TwoTempPlasmaRate +/** + * The data container `TwoTempPlasmaData` holds precalculated data common to + * all `TwoTempPlasmaRate` objects. + */ +struct TwoTempPlasmaData : public ReactionData +{ + TwoTempPlasmaData() : electronTemp(1.), logTe(0.), recipTe(1.) {} + + virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; + virtual void update(double T) override; + virtual void update(double T, double Te) override; + using ReactionData::update; + + virtual void updateTe(double Te); + + virtual void invalidateCache() override { + ReactionData::invalidateCache(); + electronTemp = NAN; + } + + double electronTemp; //!< electron temperature + double logTe; //!< logarithm of electron temperature + double recipTe; //!< inverse of electron temperature +}; + + +//! Two temperature plasma reaction rate type depends on both +//! gas temperature and electron temperature. +/*! + * The form of the two temperature plasma reaction rate coefficient is similar to an + * Arrhenius reaction rate coefficient. The temperature exponent (b) is applied to + * the electron temperature instead. In addition, the exponential term with + * activation energy for electron is included. + * + * \f[ + * k_f = A T_e^b \exp (-E_a/RT) \exp (-E_{a,e}/RT_e) + * \f] + * + * Ref.: Kossyi, I. A., Kostinsky, A. Y., Matveyev, A. A., & Silakov, V. P. (1992). + * Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures. + * Plasma Sources Science and Technology, 1(3), 207. + * doi: 10.1088/0963-0252/1/3/011 + * + * @ingroup arrheniusGroup + */ +class TwoTempPlasma : public ArrheniusBase +{ +public: + TwoTempPlasma(); + + //! Constructor. + /*! + * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units + * depend on the reaction order and the dimensionality (surface or bulk). + * @param b Temperature exponent (non-dimensional) + * @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); + + virtual const std::string type() const override { + return "two-temperature-plasma"; + } + + virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; + + //! Evaluate reaction rate + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const TwoTempPlasmaData& shared_data) const { + // m_E4_R is the electron activation (in temperature units) + return m_A * std::exp(m_b * shared_data.logTe - + m_Ea_R * shared_data.recipT + + m_E4_R * (shared_data.electronTemp - shared_data.temperature) + * shared_data.recipTe * shared_data.recipT); + } + + //! Evaluate derivative of reaction rate with respect to temperature + //! divided by reaction rate + /*! + * This method does not consider changes of electron temperature. + * A corresponding warning is raised. + * @param shared_data data shared by all reactions of a given type + */ + double ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const; + + //! Return the electron activation energy *Ea* [J/kmol] + double activationElectronEnergy() const { + return m_E4_R * GasConstant; + } +}; + +typedef BulkRate TwoTempPlasmaRate; + +} + +#endif diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 51152eccbb..e51aafe36e 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -440,20 +440,11 @@ cdef extern from "cantera/thermo/SurfPhase.h": void setCoveragesNoNorm(double*) except +translate_exception void getCoverages(double*) except +translate_exception - -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/ReactionRateFactory.h" namespace "Cantera": cdef shared_ptr[CxxReactionRate] CxxNewReactionRate "newReactionRate" (string) except +translate_exception cdef shared_ptr[CxxReactionRate] CxxNewReactionRate "newReactionRate" (CxxAnyMap&) 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 - +cdef extern from "cantera/kinetics/ReactionRate.h" namespace "Cantera": cdef cppclass CxxReactionRate "Cantera::ReactionRate": CxxReactionRate() string type() @@ -462,6 +453,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": double eval(double, vector[double]&) except +translate_exception CxxAnyMap parameters() except +translate_exception +cdef extern from "cantera/kinetics/Arrhenius.h" namespace "Cantera": cdef cppclass CxxArrheniusBase "Cantera::ArrheniusBase" (CxxReactionRate): CxxArrheniusBase() double preExponentialFactor() @@ -481,11 +473,6 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": CxxArrheniusRate(CxxAnyMap) except +translate_exception CxxArrheniusRate(double, double, double) - cdef cppclass CxxTwoTempPlasmaRate "Cantera::TwoTempPlasmaRate" (CxxArrheniusBase): - CxxTwoTempPlasmaRate(CxxAnyMap) except +translate_exception - CxxTwoTempPlasmaRate(double, double, double, double) - double activationElectronEnergy() - cdef cppclass CxxBlowersMasel "Cantera::BlowersMasel" (CxxArrheniusBase): CxxBlowersMasel(double, double, double, double) double evalRate(double, double) @@ -497,6 +484,19 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": CxxBlowersMaselRate(CxxAnyMap) except +translate_exception CxxBlowersMaselRate(double, double, double, double) +cdef extern from "cantera/kinetics/TwoTempPlasmaRate.h" namespace "Cantera": + cdef cppclass CxxTwoTempPlasmaRate "Cantera::TwoTempPlasmaRate" (CxxArrheniusBase): + CxxTwoTempPlasmaRate(CxxAnyMap) except +translate_exception + CxxTwoTempPlasmaRate(double, double, double, double) + double activationElectronEnergy() + +cdef extern from "cantera/kinetics/Reaction.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 vector[shared_ptr[CxxReaction]] CxxGetReactions "getReactions" (XML_Node&) except +translate_exception + cdef vector[shared_ptr[CxxReaction]] CxxGetReactions "getReactions" (CxxAnyValue&, CxxKinetics&) except +translate_exception + cdef cppclass CxxFalloffRate "Cantera::FalloffRate" (CxxReactionRate): CxxFalloffRate() CxxFalloffRate(CxxAnyMap) except +translate_exception diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index 469bba1171..ebb6c55de4 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -4,7 +4,7 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/kinetics/Arrhenius.h" -#include "cantera/kinetics/Kinetics.h" +#include "cantera/thermo/ThermoPhase.h" namespace Cantera { @@ -126,72 +126,14 @@ void ArrheniusBase::check(const std::string& equation, const AnyMap& node) } } -TwoTempPlasma::TwoTempPlasma() - : ArrheniusBase() +bool ArrheniusData::update(const ThermoPhase& phase, const Kinetics& kin) { - m_Ea_str = "Ea-gas"; - m_E4_str = "Ea-electron"; -} - -TwoTempPlasma::TwoTempPlasma(double A, double b, double Ea, double EE) - : ArrheniusBase(A, b, Ea) -{ - m_Ea_str = "Ea-gas"; - m_E4_str = "Ea-electron"; - m_E4_R = EE / GasConstant; -} - -double TwoTempPlasma::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const -{ - warn_user("TwoTempPlasma::ddTScaledFromStruct", - "Temperature derivative does not consider changes of electron temperature."); - return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; -} - -void TwoTempPlasma::setContext(const Reaction& rxn, const Kinetics& kin) -{ - // TwoTempPlasmaReaction is for a non-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, - "TwoTempPlasmaRate does not support reversible reactions"); - } -} - -BlowersMasel::BlowersMasel() - : m_deltaH_R(0.) -{ - m_Ea_str = "Ea0"; - m_E4_str = "w"; -} - -BlowersMasel::BlowersMasel(double A, double b, double Ea0, double w) - : ArrheniusBase(A, b, Ea0) - , m_deltaH_R(0.) -{ - m_Ea_str = "Ea0"; - m_E4_str = "w"; - m_E4_R = w / GasConstant; -} - -double BlowersMasel::ddTScaledFromStruct(const BlowersMaselData& shared_data) const -{ - warn_user("BlowersMasel::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) -{ - m_stoich_coeffs.clear(); - for (const auto& sp : rxn.reactants) { - m_stoich_coeffs.emplace_back(kin.kineticsSpeciesIndex(sp.first), -sp.second); - } - for (const auto& sp : rxn.products) { - m_stoich_coeffs.emplace_back(kin.kineticsSpeciesIndex(sp.first), sp.second); + double T = phase.temperature(); + if (T == temperature) { + return false; } + update(T); + return true; } } diff --git a/src/kinetics/BlowersMaselRate.cpp b/src/kinetics/BlowersMaselRate.cpp new file mode 100644 index 0000000000..81ed5f2211 --- /dev/null +++ b/src/kinetics/BlowersMaselRate.cpp @@ -0,0 +1,81 @@ +//! @file BlowersMaselRate.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/BlowersMaselRate.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/Kinetics.h" +#include "cantera/thermo/ThermoPhase.h" + +namespace Cantera +{ + +BlowersMaselData::BlowersMaselData() + : ready(false) + , density(NAN) + , m_state_mf_number(-1) +{ +} + +void BlowersMaselData::update(double T) { + warn_user("BlowersMaselData::update", + "This method does not update the change of reaction enthalpy."); + ReactionData::update(T); +} + +bool BlowersMaselData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + double rho = phase.density(); + int mf = phase.stateMFNumber(); + double T = phase.temperature(); + bool changed = false; + if (T != temperature) { + ReactionData::update(T); + changed = true; + } + if (changed || rho != density || mf != m_state_mf_number) { + density = rho; + m_state_mf_number = mf; + phase.getPartialMolarEnthalpies(partialMolarEnthalpies.data()); + changed = true; + } + return changed; +} + +BlowersMasel::BlowersMasel() + : m_deltaH_R(0.) +{ + m_Ea_str = "Ea0"; + m_E4_str = "w"; +} + +BlowersMasel::BlowersMasel(double A, double b, double Ea0, double w) + : ArrheniusBase(A, b, Ea0) + , m_deltaH_R(0.) +{ + m_Ea_str = "Ea0"; + m_E4_str = "w"; + m_E4_R = w / GasConstant; +} + +double BlowersMasel::ddTScaledFromStruct(const BlowersMaselData& shared_data) const +{ + warn_user("BlowersMasel::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) +{ + m_stoich_coeffs.clear(); + for (const auto& sp : rxn.reactants) { + m_stoich_coeffs.emplace_back(kin.kineticsSpeciesIndex(sp.first), -sp.second); + } + for (const auto& sp : rxn.products) { + m_stoich_coeffs.emplace_back(kin.kineticsSpeciesIndex(sp.first), sp.second); + } +} + +} diff --git a/src/kinetics/ChebyshevRate.cpp b/src/kinetics/ChebyshevRate.cpp new file mode 100644 index 0000000000..dc7b94bd17 --- /dev/null +++ b/src/kinetics/ChebyshevRate.cpp @@ -0,0 +1,174 @@ +//! @file ChebyshevRate.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/ChebyshevRate.h" +#include "cantera/thermo/ThermoPhase.h" + +namespace Cantera +{ + +void ChebyshevData::update(double T) +{ + throw CanteraError("ChebyshevData::update", + "Missing state information: 'ChebyshevData' requires pressure."); +} + +bool ChebyshevData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + double T = phase.temperature(); + double P = phase.pressure(); + if (P != pressure || T != temperature) { + update(T, P); + return true; + } + return false; +} + +void ChebyshevData::perturbPressure(double deltaP) +{ + if (m_pressure_buf > 0.) { + throw CanteraError("ChebyshevData::perturbPressure", + "Cannot apply another perturbation as state is already perturbed."); + } + m_pressure_buf = pressure; + update(temperature, pressure * (1. + deltaP)); +} + +void ChebyshevData::restore() +{ + ReactionData::restore(); + // only restore if there is a valid buffered value + if (m_pressure_buf < 0.) { + return; + } + update(temperature, m_pressure_buf); + m_pressure_buf = -1.; +} + +ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs) : ChebyshevRate() +{ + setLimits(Tmin, Tmax, Pmin, Pmax); + setData(coeffs); +} + +void ChebyshevRate::setParameters(const AnyMap& node, const UnitStack& units) +{ + m_rate_units = units.product(); + const UnitSystem& unit_system = node.units(); + Array2D coeffs; + if (node.hasKey("data")) { + const auto& T_range = node["temperature-range"].asVector(2); + const auto& P_range = node["pressure-range"].asVector(2); + auto& vcoeffs = node["data"].asVector(); + coeffs = Array2D(vcoeffs.size(), vcoeffs[0].size()); + for (size_t i = 0; i < coeffs.nRows(); i++) { + if (vcoeffs[i].size() != vcoeffs[0].size()) { + throw InputFileError("ChebyshevRate::setParameters", node["data"], + "Inconsistent number of coefficients in row {} of matrix", i + 1); + } + for (size_t j = 0; j < coeffs.nColumns(); j++) { + coeffs(i, j) = vcoeffs[i][j]; + } + } + if (m_rate_units.factor()) { + coeffs(0, 0) += std::log10(unit_system.convertTo(1.0, m_rate_units)); + } + setLimits( + unit_system.convert(T_range[0], "K"), + unit_system.convert(T_range[1], "K"), + unit_system.convert(P_range[0], "Pa"), + unit_system.convert(P_range[1], "Pa") + ); + } else { + // ensure that reaction rate can be evaluated (but returns NaN) + coeffs = Array2D(1, 1); + coeffs(0, 0) = NAN; + setLimits(290., 3000., 1.e-7, 1.e14); + } + + setData(coeffs); +} + +void ChebyshevRate::setup(double Tmin, double Tmax, double Pmin, double Pmax, + const Array2D& coeffs) +{ + warn_deprecated("ChebyshevRate::setup", "Deprecated in Cantera 2.6; " + "replaceable with setLimits() and setData()."); + setLimits(Tmin, Tmax, Pmin, Pmax); + setData(coeffs); +} + +void ChebyshevRate::setLimits(double Tmin, double Tmax, double Pmin, double Pmax) +{ + double logPmin = std::log10(Pmin); + double logPmax = std::log10(Pmax); + double TminInv = 1.0 / Tmin; + double TmaxInv = 1.0 / Tmax; + + TrNum_ = - TminInv - TmaxInv; + TrDen_ = 1.0 / (TmaxInv - TminInv); + PrNum_ = - logPmin - logPmax; + PrDen_ = 1.0 / (logPmax - logPmin); + + Tmin_ = Tmin; + Tmax_ = Tmax; + Pmin_ = Pmin; + Pmax_ = Pmax; +} + +void ChebyshevRate::setData(const Array2D& coeffs) +{ + m_coeffs = coeffs; + dotProd_.resize(coeffs.nRows()); + + // convert to row major for legacy output + // note: chebCoeffs_ is not used internally (@todo: remove after Cantera 2.6) + size_t rows = m_coeffs.nRows(); + size_t cols = m_coeffs.nColumns(); + chebCoeffs_.resize(rows * cols); + for (size_t i = 0; i < rows; i++) { + for (size_t j = 0; j < cols; j++) { + chebCoeffs_[cols * i + j] = m_coeffs(i, j); + } + } +} + +void ChebyshevRate::getParameters(AnyMap& rateNode) const +{ + rateNode["type"] = type(); + if (!m_coeffs.data().size() || std::isnan(m_coeffs(0, 0))) { + // object not fully set up + return; + } + rateNode["temperature-range"].setQuantity({Tmin(), Tmax()}, "K"); + rateNode["pressure-range"].setQuantity({Pmin(), Pmax()}, "Pa"); + size_t nT = m_coeffs.nRows(); + size_t nP = m_coeffs.nColumns(); + std::vector coeffs2d(nT, vector_fp(nP)); + for (size_t i = 0; i < nT; i++) { + for (size_t j = 0; j < nP; j++) { + coeffs2d[i][j] = m_coeffs(i, j); + } + } + // Unit conversions must take place later, after the destination unit system + // is known. A lambda function is used here to override the default behavior + Units rate_units2 = m_rate_units; + auto converter = [rate_units2](AnyValue& coeffs, const UnitSystem& units) { + if (rate_units2.factor() != 0.0) { + coeffs.asVector()[0][0] += \ + std::log10(units.convertFrom(1.0, rate_units2)); + } else if (units.getDelta(UnitSystem()).size()) { + throw CanteraError("ChebyshevRate::getParameters lambda", + "Cannot convert rate constant with unknown dimensions to a " + "non-default unit system"); + } + }; + AnyValue coeffs; + coeffs = std::move(coeffs2d); + rateNode["data"].setQuantity(coeffs, converter); +} + +} diff --git a/src/kinetics/Falloff.cpp b/src/kinetics/Falloff.cpp index e06e5862b9..414da9232b 100644 --- a/src/kinetics/Falloff.cpp +++ b/src/kinetics/Falloff.cpp @@ -11,11 +11,77 @@ #include "cantera/base/global.h" #include "cantera/base/AnyMap.h" #include "cantera/kinetics/Falloff.h" - +#include "cantera/kinetics/Kinetics.h" +#include "cantera/thermo/ThermoPhase.h" namespace Cantera { +FalloffData::FalloffData() + : ready(false) + , molar_density(NAN) + , m_state_mf_number(-1) + , m_perturbed(false) +{ + conc_3b.resize(1, NAN); + m_conc_3b_buf.resize(1, NAN); +} + +void FalloffData::update(double T) +{ + throw CanteraError("FalloffData::update", + "Missing state information: 'FalloffData' requires third-body concentration."); +} + +void FalloffData::update(double T, double M) +{ + ReactionData::update(T); + conc_3b[0] = M; +} + +bool FalloffData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + double rho_m = phase.molarDensity(); + int mf = phase.stateMFNumber(); + double T = phase.temperature(); + bool changed = false; + if (T != temperature) { + ReactionData::update(T); + changed = true; + } + if (rho_m != molar_density || mf != m_state_mf_number) { + molar_density = rho_m; + m_state_mf_number = mf; + conc_3b = kin.thirdBodyConcentrations(); + changed = true; + } + return changed; +} + +void FalloffData::perturbThirdBodies(double deltaM) +{ + if (m_perturbed) { + throw CanteraError("FalloffData::perturbThirdBodies", + "Cannot apply another perturbation as state is already perturbed."); + } + m_conc_3b_buf = conc_3b; + for (auto& c3b : conc_3b) { + c3b *= 1. + deltaM; + } + m_perturbed = true; +} + +void FalloffData::restore() +{ + ReactionData::restore(); + // only restore if there is a valid buffered value + if (!m_perturbed) { + return; + } + conc_3b = m_conc_3b_buf; + m_perturbed = false; +} + void FalloffRate::init(const vector_fp& c) { setFalloffCoeffs(c); diff --git a/src/kinetics/InterfaceRate.cpp b/src/kinetics/InterfaceRate.cpp index 48bfcd4386..842c8dc8c4 100644 --- a/src/kinetics/InterfaceRate.cpp +++ b/src/kinetics/InterfaceRate.cpp @@ -5,12 +5,95 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/Kinetics.h" +#include "cantera/kinetics/Reaction.h" #include "cantera/thermo/ThermoPhase.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/base/AnyMap.h" namespace Cantera { + +InterfaceData::InterfaceData() + : sqrtT(NAN) +{ +} + +void InterfaceData::update(double T) +{ + throw CanteraError("InterfaceData::update", + "Missing state information: 'InterfaceData' requires species coverages."); +} + +void InterfaceData::update(double T, const vector_fp& values) +{ + warn_user("InterfaceData::update", + "This method does not update the site density."); + ReactionData::update(T); + sqrtT = sqrt(T); + if (coverages.size() == 0) { + coverages = values; + logCoverages.resize(values.size()); + } else if (values.size() == coverages.size()) { + std::copy(values.begin(), values.end(), coverages.begin()); + } else { + throw CanteraError("InterfaceData::update", + "Incompatible lengths of coverage arrays: received {} elements while " + "{} are required.", values.size(), coverages.size()); + } + for (size_t n = 0; n < coverages.size(); n++) { + logCoverages[n] = std::log(std::max(coverages[n], Tiny)); + } +} + +bool InterfaceData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + int mf = 0; + for (size_t n = 0; n < kin.nPhases(); n++) { + mf += kin.thermo(n).stateMFNumber(); + } + + double T = phase.temperature(); + bool changed = false; + const auto& surf = dynamic_cast( + kin.thermo(kin.surfacePhaseIndex())); + double site_density = surf.siteDensity(); + if (density != site_density) { + density = surf.siteDensity(); + changed = true; + } + if (T != temperature) { + ReactionData::update(T); + sqrtT = sqrt(T); + changed = true; + } + if (changed || mf != m_state_mf_number) { + surf.getCoverages(coverages.data()); + for (size_t n = 0; n < coverages.size(); n++) { + logCoverages[n] = std::log(std::max(coverages[n], Tiny)); + } + for (size_t n = 0; n < kin.nPhases(); n++) { + size_t start = kin.kineticsSpeciesIndex(0, n); + const auto& ph = kin.thermo(n); + electricPotentials[n] = ph.electricPotential(); + ph.getPartialMolarEnthalpies(partialMolarEnthalpies.data() + start); + ph.getStandardChemPotentials(standardChemPotentials.data() + start); + size_t nsp = ph.nSpecies(); + for (size_t k = 0; k < nsp; k++) { + // only used for exchange current density formulation + standardConcentrations[k + start] = ph.standardConcentration(k); + } + } + m_state_mf_number = mf; + changed = true; + } + return changed; +} + +void InterfaceData::perturbTemperature(double deltaT) +{ + throw NotImplementedError("InterfaceData::perturbTemperature"); +} + InterfaceRateBase::InterfaceRateBase() : m_siteDensity(NAN) , m_acov(0.) diff --git a/src/kinetics/PlogRate.cpp b/src/kinetics/PlogRate.cpp new file mode 100644 index 0000000000..44c5a51e56 --- /dev/null +++ b/src/kinetics/PlogRate.cpp @@ -0,0 +1,208 @@ +//! @file PlogRate.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/PlogRate.h" +#include "cantera/kinetics/RxnRates.h" +#include "cantera/thermo/ThermoPhase.h" + +namespace Cantera +{ + +void PlogData::update(double T) +{ + throw CanteraError("PlogData::update", + "Missing state information: 'PlogData' requires pressure."); +} + +bool PlogData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + double T = phase.temperature(); + double P = phase.pressure(); + if (P != pressure || T != temperature) { + update(T, P); + return true; + } + return false; +} + +void PlogData::perturbPressure(double deltaP) +{ + if (m_pressure_buf > 0.) { + throw CanteraError("PlogData::perturbPressure", + "Cannot apply another perturbation as state is already perturbed."); + } + m_pressure_buf = pressure; + update(temperature, pressure * (1. + deltaP)); +} + +void PlogData::restore() +{ + ReactionData::restore(); + // only restore if there is a valid buffered value + if (m_pressure_buf < 0.) { + return; + } + update(temperature, m_pressure_buf); + m_pressure_buf = -1.; +} + +PlogRate::PlogRate() + : logP_(-1000) + , logP1_(1000) + , logP2_(-1000) + , rDeltaP_(-1.0) +{ +} + +PlogRate::PlogRate(const std::multimap& rates) + : PlogRate() +{ + setRates(rates); +} + +PlogRate::PlogRate(const std::multimap& rates) + : PlogRate() +{ + setup(rates); +} + +void PlogRate::setParameters(const AnyMap& node, const UnitStack& units) +{ + auto rate_units = units.product(); + 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()) { + for (const auto& rate : rates) { + multi_rates.insert({rate.convert("P", "Pa"), + Arrhenius3(AnyValue(rate), units, rate_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); +} + +void PlogRate::getParameters(AnyMap& rateNode, const Units& rate_units) const +{ + std::vector rateList; + rateNode["type"] = type(); + if (!rates_.size() + || (rates_.size() > 1 && std::isnan(rates_[1].preExponentialFactor()))) + { + // object not fully set up + return; + } + for (const auto& r : getRates()) { + AnyMap rateNode_; + rateNode_["P"].setQuantity(r.first, "Pa"); + if (rate_units.factor() == 0) { + r.second.getRateParameters(rateNode_); + } else { + // legacy rate + Arrhenius2(r.second).getParameters(rateNode_, rate_units); + } + rateList.push_back(std::move(rateNode_)); + } + rateNode["rate-constants"] = std::move(rateList); +} + +void PlogRate::setup(const std::multimap& rates) +{ + std::multimap rates2; + for (const auto& item : rates) { + rates2.emplace(item.first, item.second); + } + setRates(rates2); +} + +void PlogRate::setRates(const std::multimap& rates) +{ + size_t j = 0; + rates_.clear(); + pressures_.clear(); + rates_.reserve(rates.size()); + // Insert intermediate pressures + for (const auto& rate : rates) { + double logp = std::log(rate.first); + if (pressures_.empty() || pressures_.rbegin()->first != logp) { + // starting a new group + pressures_[logp] = {j, j+1}; + } else { + // another rate expression at the same pressure + pressures_[logp].second = j+1; + } + + j++; + rates_.push_back(rate.second); + } + + // Duplicate the first and last groups to handle P < P_0 and P > P_N + pressures_.insert({-1000.0, pressures_.begin()->second}); + pressures_.insert({1000.0, pressures_.rbegin()->second}); +} + +void PlogRate::validate(const std::string& equation) +{ + fmt::memory_buffer err_reactions; + double T[] = {200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0}; + + for (auto iter = ++pressures_.begin(); iter->first < 1000; iter++) { + update_C(&iter->first); + for (size_t i=0; i < 6; i++) { + double k = 0; + for (size_t p = ilow1_; p < ilow2_; p++) { + k += rates_.at(p).evalRate(log(T[i]), 1.0 / T[i]); + } + if (k < 0) { + fmt_append(err_reactions, + "\nInvalid rate coefficient for reaction '{}'\n" + "at P = {:.5g}, T = {:.1f}\n", + equation, std::exp(iter->first), T[i]); + } + } + } + if (err_reactions.size()) { + throw CanteraError("PlogRate::validate", to_string(err_reactions)); + } +} + +std::vector > PlogRate::rates() const +{ + auto rateMap = getRates(); + std::vector> out; + for (const auto& item : rateMap) { + out.emplace_back(item.first, Arrhenius2(item.second)); + } + return out; +} + +std::multimap PlogRate::getRates() const +{ + std::multimap rateMap; + // initial preincrement to skip rate for P --> 0 + for (auto iter = ++pressures_.begin(); + iter->first < 1000; // skip rates for (P --> infinity) + ++iter) { + for (size_t i = iter->second.first; i < iter->second.second; i++) { + rateMap.insert({std::exp(iter->first), rates_[i]}); + } + } + return rateMap; +} + +} diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index ee90e680c6..0f9188d288 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -6,7 +6,7 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/kinetics/Reaction.h" -#include "cantera/kinetics/ReactionFactory.h" +#include "ReactionFactory.h" #include "cantera/kinetics/ReactionRateFactory.h" #include "cantera/kinetics/FalloffFactory.h" #include "cantera/kinetics/Kinetics.h" @@ -1261,6 +1261,105 @@ CustomFunc1Reaction::CustomFunc1Reaction(const AnyMap& node, const Kinetics& kin } } +bool isThreeBody(const Reaction& R) +{ + // detect explicitly specified collision partner + size_t found = 0; + for (const auto& reac : R.reactants) { + auto prod = R.products.find(reac.first); + if (prod != R.products.end() && + trunc(reac.second) == reac.second && trunc(prod->second) == prod->second) { + // candidate species with integer stoichiometric coefficients on both sides + found += 1; + } + } + if (found != 1) { + return false; + } + + // ensure that all reactants have integer stoichiometric coefficients + size_t nreac = 0; + for (const auto& reac : R.reactants) { + if (trunc(reac.second) != reac.second) { + return false; + } + nreac += static_cast(reac.second); + } + + // ensure that all products have integer stoichiometric coefficients + size_t nprod = 0; + for (const auto& prod : R.products) { + if (trunc(prod.second) != prod.second) { + return false; + } + nprod += static_cast(prod.second); + } + + // either reactant or product side involves exactly three species + return (nreac == 3) || (nprod == 3); +} + +unique_ptr newReaction(const std::string& type) +{ + AnyMap rxn_node; + Kinetics kin; + unique_ptr R(ReactionFactory::factory()->create(type, rxn_node, kin)); + return R; +} + +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"; + } + + if (!(ReactionFactoryXML::factory()->exists(type))) { + throw CanteraError("newReaction", + "Unknown reaction type '" + rxn_node["type"] + "'"); + } + if (type.empty()) { + // Reaction type is not specified + // See if this is a three-body reaction with a specified collision partner + ElementaryReaction2 testReaction; + setupReaction(testReaction, rxn_node); + if (isThreeBody(testReaction)) { + type = "three-body"; + } + } + Reaction* R = ReactionFactoryXML::factory()->create(type, rxn_node); + return unique_ptr(R); +} + +unique_ptr newReaction(const AnyMap& rxn_node, const Kinetics& kin) +{ + std::string type = "elementary"; + size_t nDim = kin.thermo(kin.reactionPhaseIndex()).nDim(); + if (rxn_node.hasKey("type")) { + type = rxn_node["type"].asString(); + } else if (nDim == 3) { + // Reaction type is not specified + // See if this is a three-body reaction with a specified collision partner + ElementaryReaction2 testReaction; + parseReactionEquation(testReaction, rxn_node["equation"].asString(), + rxn_node, &kin); + if (isThreeBody(testReaction)) { + type = "three-body"; + } + } + + if (!(ReactionFactory::factory()->exists(type))) { + throw InputFileError("ReactionFactory::newReaction", rxn_node["type"], + "Unknown reaction type '{}'", type); + } + Reaction* R = ReactionFactory::factory()->create(type, rxn_node, kin); + return unique_ptr(R); +} + Arrhenius2 readArrhenius(const XML_Node& arrhenius_node) { return Arrhenius2(getFloat(arrhenius_node, "A", "toSI"), diff --git a/src/kinetics/ReactionData.cpp b/src/kinetics/ReactionData.cpp deleted file mode 100644 index be18f4387f..0000000000 --- a/src/kinetics/ReactionData.cpp +++ /dev/null @@ -1,346 +0,0 @@ -//! @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/kinetics/Kinetics.h" -#include "cantera/thermo/ThermoPhase.h" -#include "cantera/thermo/SurfPhase.h" -#include "cantera/base/ctexceptions.h" - -namespace Cantera -{ - -void ReactionData::update(double T, double extra) -{ - throw NotImplementedError("ReactionData::update", - "ReactionData type does not use extra scalar argument."); -} - -void ReactionData::update(double T, const vector_fp& extra) -{ - throw NotImplementedError("ReactionData::update", - "ReactionData type does not use extra vector argument."); -} - -void ReactionData::perturbTemperature(double deltaT) -{ - if (m_temperature_buf > 0.) { - throw CanteraError("ReactionData::perturbTemperature", - "Cannot apply another perturbation as state is already perturbed."); - } - m_temperature_buf = temperature; - ReactionData::update(temperature * (1. + deltaT)); -} - -void ReactionData::restore() -{ - // only restore if there is a valid buffered value - if (m_temperature_buf < 0.) { - return; - } - ReactionData::update(m_temperature_buf); - m_temperature_buf = -1.; -} - -bool ArrheniusData::update(const ThermoPhase& phase, const Kinetics& kin) -{ - double T = phase.temperature(); - if (T == temperature) { - return false; - } - update(T); - return true; -} - -bool TwoTempPlasmaData::update(const ThermoPhase& phase, const Kinetics& kin) -{ - double T = phase.temperature(); - double Te = phase.electronTemperature(); - bool changed = false; - if (T != temperature) { - ReactionData::update(T); - changed = true; - } - if (Te != electronTemp) { - updateTe(Te); - changed = true; - } - return changed; -} - -void TwoTempPlasmaData::update(double T) -{ - throw CanteraError("TwoTempPlasmaData::update", - "Missing state information: 'TwoTempPlasmaData' requires electron temperature."); -} - -void TwoTempPlasmaData::update(double T, double Te) -{ - ReactionData::update(T); - updateTe(Te); -} - -void TwoTempPlasmaData::updateTe(double Te) -{ - electronTemp = Te; - logTe = std::log(Te); - recipTe = 1./Te; -} - -BlowersMaselData::BlowersMaselData() - : ready(false) - , density(NAN) - , m_state_mf_number(-1) -{ -} - -void BlowersMaselData::update(double T) { - warn_user("BlowersMaselData::update", - "This method does not update the change of reaction enthalpy."); - ReactionData::update(T); -} - -bool BlowersMaselData::update(const ThermoPhase& phase, const Kinetics& kin) -{ - double rho = phase.density(); - int mf = phase.stateMFNumber(); - double T = phase.temperature(); - bool changed = false; - if (T != temperature) { - ReactionData::update(T); - changed = true; - } - if (changed || rho != density || mf != m_state_mf_number) { - density = rho; - m_state_mf_number = mf; - phase.getPartialMolarEnthalpies(partialMolarEnthalpies.data()); - changed = true; - } - return changed; -} - -FalloffData::FalloffData() - : ready(false) - , molar_density(NAN) - , m_state_mf_number(-1) - , m_perturbed(false) -{ - conc_3b.resize(1, NAN); - m_conc_3b_buf.resize(1, NAN); -} - -void FalloffData::update(double T) -{ - throw CanteraError("FalloffData::update", - "Missing state information: 'FalloffData' requires third-body concentration."); -} - -void FalloffData::update(double T, double M) -{ - ReactionData::update(T); - conc_3b[0] = M; -} - -bool FalloffData::update(const ThermoPhase& phase, const Kinetics& kin) -{ - double rho_m = phase.molarDensity(); - int mf = phase.stateMFNumber(); - double T = phase.temperature(); - bool changed = false; - if (T != temperature) { - ReactionData::update(T); - changed = true; - } - if (rho_m != molar_density || mf != m_state_mf_number) { - molar_density = rho_m; - m_state_mf_number = mf; - conc_3b = kin.thirdBodyConcentrations(); - changed = true; - } - return changed; -} - -void FalloffData::perturbThirdBodies(double deltaM) -{ - if (m_perturbed) { - throw CanteraError("FalloffData::perturbThirdBodies", - "Cannot apply another perturbation as state is already perturbed."); - } - m_conc_3b_buf = conc_3b; - for (auto& c3b : conc_3b) { - c3b *= 1. + deltaM; - } - m_perturbed = true; -} - -void FalloffData::restore() -{ - ReactionData::restore(); - // only restore if there is a valid buffered value - if (!m_perturbed) { - return; - } - conc_3b = m_conc_3b_buf; - m_perturbed = false; -} - -void PlogData::update(double T) -{ - throw CanteraError("PlogData::update", - "Missing state information: 'PlogData' requires pressure."); -} - -bool PlogData::update(const ThermoPhase& phase, const Kinetics& kin) -{ - double T = phase.temperature(); - double P = phase.pressure(); - if (P != pressure || T != temperature) { - update(T, P); - return true; - } - return false; -} - -void PlogData::perturbPressure(double deltaP) -{ - if (m_pressure_buf > 0.) { - throw CanteraError("PlogData::perturbPressure", - "Cannot apply another perturbation as state is already perturbed."); - } - m_pressure_buf = pressure; - update(temperature, pressure * (1. + deltaP)); -} - -void PlogData::restore() -{ - ReactionData::restore(); - // only restore if there is a valid buffered value - if (m_pressure_buf < 0.) { - return; - } - update(temperature, m_pressure_buf); - m_pressure_buf = -1.; -} - -void ChebyshevData::update(double T) -{ - throw CanteraError("ChebyshevData::update", - "Missing state information: 'ChebyshevData' requires pressure."); -} - -bool ChebyshevData::update(const ThermoPhase& phase, const Kinetics& kin) -{ - double T = phase.temperature(); - double P = phase.pressure(); - if (P != pressure || T != temperature) { - update(T, P); - return true; - } - return false; -} - -void ChebyshevData::perturbPressure(double deltaP) -{ - if (m_pressure_buf > 0.) { - throw CanteraError("ChebyshevData::perturbPressure", - "Cannot apply another perturbation as state is already perturbed."); - } - m_pressure_buf = pressure; - update(temperature, pressure * (1. + deltaP)); -} - -void ChebyshevData::restore() -{ - ReactionData::restore(); - // only restore if there is a valid buffered value - if (m_pressure_buf < 0.) { - return; - } - update(temperature, m_pressure_buf); - m_pressure_buf = -1.; -} - -InterfaceData::InterfaceData() - : sqrtT(NAN) -{ -} - -void InterfaceData::update(double T) -{ - throw CanteraError("InterfaceData::update", - "Missing state information: 'InterfaceData' requires species coverages."); -} - -void InterfaceData::update(double T, const vector_fp& values) -{ - warn_user("InterfaceData::update", - "This method does not update the site density."); - ReactionData::update(T); - sqrtT = sqrt(T); - if (coverages.size() == 0) { - coverages = values; - logCoverages.resize(values.size()); - } else if (values.size() == coverages.size()) { - std::copy(values.begin(), values.end(), coverages.begin()); - } else { - throw CanteraError("InterfaceData::update", - "Incompatible lengths of coverage arrays: received {} elements while " - "{} are required.", values.size(), coverages.size()); - } - for (size_t n = 0; n < coverages.size(); n++) { - logCoverages[n] = std::log(std::max(coverages[n], Tiny)); - } -} - -bool InterfaceData::update(const ThermoPhase& phase, const Kinetics& kin) -{ - int mf = 0; - for (size_t n = 0; n < kin.nPhases(); n++) { - mf += kin.thermo(n).stateMFNumber(); - } - - double T = phase.temperature(); - bool changed = false; - const auto& surf = dynamic_cast( - kin.thermo(kin.surfacePhaseIndex())); - double site_density = surf.siteDensity(); - if (density != site_density) { - density = surf.siteDensity(); - changed = true; - } - if (T != temperature) { - ReactionData::update(T); - sqrtT = sqrt(T); - changed = true; - } - if (changed || mf != m_state_mf_number) { - surf.getCoverages(coverages.data()); - for (size_t n = 0; n < coverages.size(); n++) { - logCoverages[n] = std::log(std::max(coverages[n], Tiny)); - } - for (size_t n = 0; n < kin.nPhases(); n++) { - size_t start = kin.kineticsSpeciesIndex(0, n); - const auto& ph = kin.thermo(n); - electricPotentials[n] = ph.electricPotential(); - ph.getPartialMolarEnthalpies(partialMolarEnthalpies.data() + start); - ph.getStandardChemPotentials(standardChemPotentials.data() + start); - size_t nsp = ph.nSpecies(); - for (size_t k = 0; k < nsp; k++) { - // only used for exchange current density formulation - standardConcentrations[k + start] = ph.standardConcentration(k); - } - } - m_state_mf_number = mf; - changed = true; - } - return changed; -} - -void InterfaceData::perturbTemperature(double deltaT) -{ - throw NotImplementedError("InterfaceData::perturbTemperature"); -} - -} diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index 08c664f278..34e7dce1ee 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -7,7 +7,7 @@ #include "cantera/thermo/ThermoPhase.h" #include "cantera/kinetics/Reaction.h" -#include "cantera/kinetics/ReactionFactory.h" +#include "ReactionFactory.h" #include "cantera/kinetics/ReactionRateFactory.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/base/ctml.h" @@ -223,103 +223,4 @@ ReactionFactoryXML::ReactionFactoryXML() addAlias("interface-legacy", "electrochemical"); } -bool isThreeBody(const Reaction& R) -{ - // detect explicitly specified collision partner - size_t found = 0; - for (const auto& reac : R.reactants) { - auto prod = R.products.find(reac.first); - if (prod != R.products.end() && - trunc(reac.second) == reac.second && trunc(prod->second) == prod->second) { - // candidate species with integer stoichiometric coefficients on both sides - found += 1; - } - } - if (found != 1) { - return false; - } - - // ensure that all reactants have integer stoichiometric coefficients - size_t nreac = 0; - for (const auto& reac : R.reactants) { - if (trunc(reac.second) != reac.second) { - return false; - } - nreac += static_cast(reac.second); - } - - // ensure that all products have integer stoichiometric coefficients - size_t nprod = 0; - for (const auto& prod : R.products) { - if (trunc(prod.second) != prod.second) { - return false; - } - nprod += static_cast(prod.second); - } - - // either reactant or product side involves exactly three species - return (nreac == 3) || (nprod == 3); -} - -unique_ptr newReaction(const std::string& type) -{ - AnyMap rxn_node; - Kinetics kin; - unique_ptr R(ReactionFactory::factory()->create(type, rxn_node, kin)); - return R; -} - -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"; - } - - if (!(ReactionFactoryXML::factory()->exists(type))) { - throw CanteraError("newReaction", - "Unknown reaction type '" + rxn_node["type"] + "'"); - } - if (type.empty()) { - // Reaction type is not specified - // See if this is a three-body reaction with a specified collision partner - ElementaryReaction2 testReaction; - setupReaction(testReaction, rxn_node); - if (isThreeBody(testReaction)) { - type = "three-body"; - } - } - Reaction* R = ReactionFactoryXML::factory()->create(type, rxn_node); - return unique_ptr(R); -} - -unique_ptr newReaction(const AnyMap& rxn_node, const Kinetics& kin) -{ - std::string type = "elementary"; - size_t nDim = kin.thermo(kin.reactionPhaseIndex()).nDim(); - if (rxn_node.hasKey("type")) { - type = rxn_node["type"].asString(); - } else if (nDim == 3) { - // Reaction type is not specified - // See if this is a three-body reaction with a specified collision partner - ElementaryReaction2 testReaction; - parseReactionEquation(testReaction, rxn_node["equation"].asString(), - rxn_node, &kin); - if (isThreeBody(testReaction)) { - type = "three-body"; - } - } - - if (!(ReactionFactory::factory()->exists(type))) { - throw InputFileError("ReactionFactory::newReaction", rxn_node["type"], - "Unknown reaction type '{}'", type); - } - Reaction* R = ReactionFactory::factory()->create(type, rxn_node, kin); - return unique_ptr(R); -} - } diff --git a/include/cantera/kinetics/ReactionFactory.h b/src/kinetics/ReactionFactory.h similarity index 79% rename from include/cantera/kinetics/ReactionFactory.h rename to src/kinetics/ReactionFactory.h index db24d0dbe2..437ecdcf0a 100644 --- a/include/cantera/kinetics/ReactionFactory.h +++ b/src/kinetics/ReactionFactory.h @@ -93,26 +93,5 @@ class ReactionFactoryXML : public Factory static std::mutex reaction_mutex; }; - -//! Create a new empty Reaction object -/*! - * @param type string identifying type of reaction. - * @deprecated To be removed after Cantera 2.6. Only used for legacy reaction types. - */ -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/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 432e9d1bf0..4a44309e9f 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -10,6 +10,10 @@ #include "cantera/thermo/ThermoPhase.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/kinetics/Falloff.h" +#include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/ChebyshevRate.h" +#include "cantera/kinetics/Custom.h" +#include "cantera/kinetics/RxnRates.h" #include "cantera/kinetics/InterfaceRate.h" #include "cantera/base/AnyMap.h" diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 1a357423ef..151982dfb9 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -95,285 +95,4 @@ void SurfaceArrhenius::addCoverageDependence(size_t k, doublereal a, } } -PlogRate::PlogRate() - : logP_(-1000) - , logP1_(1000) - , logP2_(-1000) - , rDeltaP_(-1.0) -{ -} - -PlogRate::PlogRate(const std::multimap& rates) - : PlogRate() -{ - setRates(rates); -} - -PlogRate::PlogRate(const std::multimap& rates) - : PlogRate() -{ - setup(rates); -} - -void PlogRate::setParameters(const AnyMap& node, const UnitStack& units) -{ - auto rate_units = units.product(); - 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()) { - for (const auto& rate : rates) { - multi_rates.insert({rate.convert("P", "Pa"), - Arrhenius3(AnyValue(rate), units, rate_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); -} - -void PlogRate::getParameters(AnyMap& rateNode, const Units& rate_units) const -{ - std::vector rateList; - rateNode["type"] = type(); - if (!rates_.size() - || (rates_.size() > 1 && std::isnan(rates_[1].preExponentialFactor()))) - { - // object not fully set up - return; - } - for (const auto& r : getRates()) { - AnyMap rateNode_; - rateNode_["P"].setQuantity(r.first, "Pa"); - if (rate_units.factor() == 0) { - r.second.getRateParameters(rateNode_); - } else { - // legacy rate - Arrhenius2(r.second).getParameters(rateNode_, rate_units); - } - rateList.push_back(std::move(rateNode_)); - } - rateNode["rate-constants"] = std::move(rateList); -} - -void PlogRate::setup(const std::multimap& rates) -{ - std::multimap rates2; - for (const auto& item : rates) { - rates2.emplace(item.first, item.second); - } - setRates(rates2); -} - -void PlogRate::setRates(const std::multimap& rates) -{ - size_t j = 0; - rates_.clear(); - pressures_.clear(); - rates_.reserve(rates.size()); - // Insert intermediate pressures - for (const auto& rate : rates) { - double logp = std::log(rate.first); - if (pressures_.empty() || pressures_.rbegin()->first != logp) { - // starting a new group - pressures_[logp] = {j, j+1}; - } else { - // another rate expression at the same pressure - pressures_[logp].second = j+1; - } - - j++; - rates_.push_back(rate.second); - } - - // Duplicate the first and last groups to handle P < P_0 and P > P_N - pressures_.insert({-1000.0, pressures_.begin()->second}); - pressures_.insert({1000.0, pressures_.rbegin()->second}); -} - -void PlogRate::validate(const std::string& equation) -{ - fmt::memory_buffer err_reactions; - double T[] = {200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0}; - - for (auto iter = ++pressures_.begin(); iter->first < 1000; iter++) { - update_C(&iter->first); - for (size_t i=0; i < 6; i++) { - double k = 0; - for (size_t p = ilow1_; p < ilow2_; p++) { - k += rates_.at(p).evalRate(log(T[i]), 1.0 / T[i]); - } - if (k < 0) { - fmt_append(err_reactions, - "\nInvalid rate coefficient for reaction '{}'\n" - "at P = {:.5g}, T = {:.1f}\n", - equation, std::exp(iter->first), T[i]); - } - } - } - if (err_reactions.size()) { - throw CanteraError("PlogRate::validate", to_string(err_reactions)); - } -} - -std::vector > PlogRate::rates() const -{ - auto rateMap = getRates(); - std::vector> out; - for (const auto& item : rateMap) { - out.emplace_back(item.first, Arrhenius2(item.second)); - } - return out; -} - -std::multimap PlogRate::getRates() const -{ - std::multimap rateMap; - // initial preincrement to skip rate for P --> 0 - for (auto iter = ++pressures_.begin(); - iter->first < 1000; // skip rates for (P --> infinity) - ++iter) { - for (size_t i = iter->second.first; i < iter->second.second; i++) { - rateMap.insert({std::exp(iter->first), rates_[i]}); - } - } - return rateMap; -} - -ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, - const Array2D& coeffs) : ChebyshevRate() -{ - setLimits(Tmin, Tmax, Pmin, Pmax); - setData(coeffs); -} - -void ChebyshevRate::setParameters(const AnyMap& node, const UnitStack& units) -{ - m_rate_units = units.product(); - const UnitSystem& unit_system = node.units(); - Array2D coeffs; - if (node.hasKey("data")) { - const auto& T_range = node["temperature-range"].asVector(2); - const auto& P_range = node["pressure-range"].asVector(2); - auto& vcoeffs = node["data"].asVector(); - coeffs = Array2D(vcoeffs.size(), vcoeffs[0].size()); - for (size_t i = 0; i < coeffs.nRows(); i++) { - if (vcoeffs[i].size() != vcoeffs[0].size()) { - throw InputFileError("ChebyshevRate::setParameters", node["data"], - "Inconsistent number of coefficients in row {} of matrix", i + 1); - } - for (size_t j = 0; j < coeffs.nColumns(); j++) { - coeffs(i, j) = vcoeffs[i][j]; - } - } - if (m_rate_units.factor()) { - coeffs(0, 0) += std::log10(unit_system.convertTo(1.0, m_rate_units)); - } - setLimits( - unit_system.convert(T_range[0], "K"), - unit_system.convert(T_range[1], "K"), - unit_system.convert(P_range[0], "Pa"), - unit_system.convert(P_range[1], "Pa") - ); - } else { - // ensure that reaction rate can be evaluated (but returns NaN) - coeffs = Array2D(1, 1); - coeffs(0, 0) = NAN; - setLimits(290., 3000., 1.e-7, 1.e14); - } - - setData(coeffs); -} - -void ChebyshevRate::setup(double Tmin, double Tmax, double Pmin, double Pmax, - const Array2D& coeffs) -{ - warn_deprecated("ChebyshevRate::setup", "Deprecated in Cantera 2.6; " - "replaceable with setLimits() and setData()."); - setLimits(Tmin, Tmax, Pmin, Pmax); - setData(coeffs); -} - -void ChebyshevRate::setLimits(double Tmin, double Tmax, double Pmin, double Pmax) -{ - double logPmin = std::log10(Pmin); - double logPmax = std::log10(Pmax); - double TminInv = 1.0 / Tmin; - double TmaxInv = 1.0 / Tmax; - - TrNum_ = - TminInv - TmaxInv; - TrDen_ = 1.0 / (TmaxInv - TminInv); - PrNum_ = - logPmin - logPmax; - PrDen_ = 1.0 / (logPmax - logPmin); - - Tmin_ = Tmin; - Tmax_ = Tmax; - Pmin_ = Pmin; - Pmax_ = Pmax; -} - -void ChebyshevRate::setData(const Array2D& coeffs) -{ - m_coeffs = coeffs; - dotProd_.resize(coeffs.nRows()); - - // convert to row major for legacy output - // note: chebCoeffs_ is not used internally (@todo: remove after Cantera 2.6) - size_t rows = m_coeffs.nRows(); - size_t cols = m_coeffs.nColumns(); - chebCoeffs_.resize(rows * cols); - for (size_t i = 0; i < rows; i++) { - for (size_t j = 0; j < cols; j++) { - chebCoeffs_[cols * i + j] = m_coeffs(i, j); - } - } -} - -void ChebyshevRate::getParameters(AnyMap& rateNode) const -{ - rateNode["type"] = type(); - if (!m_coeffs.data().size() || std::isnan(m_coeffs(0, 0))) { - // object not fully set up - return; - } - rateNode["temperature-range"].setQuantity({Tmin(), Tmax()}, "K"); - rateNode["pressure-range"].setQuantity({Pmin(), Pmax()}, "Pa"); - size_t nT = m_coeffs.nRows(); - size_t nP = m_coeffs.nColumns(); - std::vector coeffs2d(nT, vector_fp(nP)); - for (size_t i = 0; i < nT; i++) { - for (size_t j = 0; j < nP; j++) { - coeffs2d[i][j] = m_coeffs(i, j); - } - } - // Unit conversions must take place later, after the destination unit system - // is known. A lambda function is used here to override the default behavior - Units rate_units2 = m_rate_units; - auto converter = [rate_units2](AnyValue& coeffs, const UnitSystem& units) { - if (rate_units2.factor() != 0.0) { - coeffs.asVector()[0][0] += \ - std::log10(units.convertFrom(1.0, rate_units2)); - } else if (units.getDelta(UnitSystem()).size()) { - throw CanteraError("ChebyshevRate::getParameters lambda", - "Cannot convert rate constant with unknown dimensions to a " - "non-default unit system"); - } - }; - AnyValue coeffs; - coeffs = std::move(coeffs2d); - rateNode["data"].setQuantity(coeffs, converter); -} - } diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp new file mode 100644 index 0000000000..7b5a2a0f37 --- /dev/null +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -0,0 +1,81 @@ +//! @file TwoTempPlasmaRate.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/TwoTempPlasmaRate.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/thermo/ThermoPhase.h" + +namespace Cantera +{ + +bool TwoTempPlasmaData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + double T = phase.temperature(); + double Te = phase.electronTemperature(); + bool changed = false; + if (T != temperature) { + ReactionData::update(T); + changed = true; + } + if (Te != electronTemp) { + updateTe(Te); + changed = true; + } + return changed; +} + +void TwoTempPlasmaData::update(double T) +{ + throw CanteraError("TwoTempPlasmaData::update", + "Missing state information: 'TwoTempPlasmaData' requires electron temperature."); +} + +void TwoTempPlasmaData::update(double T, double Te) +{ + ReactionData::update(T); + updateTe(Te); +} + +void TwoTempPlasmaData::updateTe(double Te) +{ + electronTemp = Te; + logTe = std::log(Te); + recipTe = 1./Te; +} + +TwoTempPlasma::TwoTempPlasma() + : ArrheniusBase() +{ + m_Ea_str = "Ea-gas"; + m_E4_str = "Ea-electron"; +} + +TwoTempPlasma::TwoTempPlasma(double A, double b, double Ea, double EE) + : ArrheniusBase(A, b, Ea) +{ + m_Ea_str = "Ea-gas"; + m_E4_str = "Ea-electron"; + m_E4_R = EE / GasConstant; +} + +double TwoTempPlasma::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const +{ + warn_user("TwoTempPlasma::ddTScaledFromStruct", + "Temperature derivative does not consider changes of electron temperature."); + return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; +} + +void TwoTempPlasma::setContext(const Reaction& rxn, const Kinetics& kin) +{ + // TwoTempPlasmaReaction is for a non-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, + "TwoTempPlasmaRate does not support reversible reactions"); + } +} + +} diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index 219567f0b6..05777ffd99 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -16,7 +16,6 @@ #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/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 8c1dc91950..96d37cbebf 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -8,6 +8,7 @@ #include "cantera/zeroD/Wall.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/thermo/SurfPhase.h" +#include "cantera/base/utilities.h" using namespace std; diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 8418f08ec7..16f509b86c 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -8,6 +8,7 @@ #include "cantera/zeroD/ReactorNet.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/thermo/ThermoPhase.h" +#include "cantera/base/utilities.h" using namespace std; diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 14f25419e5..a79b1b6fc7 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -8,6 +8,7 @@ #include "cantera/zeroD/Wall.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/thermo/ThermoPhase.h" +#include "cantera/base/utilities.h" using namespace std; diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 20fe4577cf..99554fa4b5 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -11,6 +11,7 @@ #include "cantera/zeroD/ReactorSurface.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/base/Solution.h" +#include "cantera/base/utilities.h" #include diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index d7923b279b..5ab3ff1f8b 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -3,9 +3,9 @@ #include "cantera/base/Solution.h" #include "cantera/base/Interface.h" #include "cantera/kinetics/GasKinetics.h" +#include "cantera/kinetics/TwoTempPlasmaRate.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/kinetics/KineticsFactory.h" -#include "cantera/kinetics/ReactionFactory.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/base/Array.h" diff --git a/test/kinetics/pdep.cpp b/test/kinetics/pdep.cpp index fd89b20680..c30a437e7d 100644 --- a/test/kinetics/pdep.cpp +++ b/test/kinetics/pdep.cpp @@ -2,6 +2,7 @@ #include "cantera/kinetics/Kinetics.h" #include "cantera/thermo/ThermoPhase.h" #include "cantera/base/Solution.h" +#include "cantera/base/global.h" namespace Cantera {