From 3b71d75ada6202cd4620e65f47396dfb75dfcfba Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 30 Oct 2014 21:09:59 +0000 Subject: [PATCH] [Kinetics] Add class BulkKinetics This class serves as a common base for GasKinetics and AqueousKinetics, eliminating most of the redundancy between the two classes. --- include/cantera/kinetics/AqueousKinetics.h | 130 +---------- include/cantera/kinetics/BulkKinetics.h | 63 +++++ include/cantera/kinetics/GasKinetics.h | 61 +---- src/kinetics/AqueousKinetics.cpp | 249 +------------------- src/kinetics/BulkKinetics.cpp | 164 +++++++++++++ src/kinetics/GasKinetics.cpp | 257 ++------------------- 6 files changed, 262 insertions(+), 662 deletions(-) create mode 100644 include/cantera/kinetics/BulkKinetics.h create mode 100644 src/kinetics/BulkKinetics.cpp diff --git a/include/cantera/kinetics/AqueousKinetics.h b/include/cantera/kinetics/AqueousKinetics.h index 1245c5ec1..9bc8dbf95 100644 --- a/include/cantera/kinetics/AqueousKinetics.h +++ b/include/cantera/kinetics/AqueousKinetics.h @@ -1,6 +1,5 @@ /** * @file AqueousKinetics.h - * * @ingroup chemkinetics */ @@ -9,44 +8,23 @@ #ifndef CT_AQUEOUSKINETICS_H #define CT_AQUEOUSKINETICS_H -#include "Kinetics.h" -#include "ReactionStoichMgr.h" -#include "RateCoeffMgr.h" -#include "cantera/base/utilities.h" +#include "BulkKinetics.h" namespace Cantera { -// forward references - -class ReactionData; - /** - * Kinetics manager for elementary aqueous-phase chemistry. This - * kinetics manager implements standard mass-action reaction rate - * expressions for liquids - * - * - * Concentration + * Kinetics manager for elementary aqueous-phase chemistry. This kinetics + * manager implements standard mass-action reaction rate expressions for liquids * * @ingroup kinetics - * @deprecated Not actually implemented */ -class AqueousKinetics : public Kinetics +class AqueousKinetics : public BulkKinetics { - public: - - //! @name Constructors - //! @{ - /// Constructor. Creates an empty reaction mechanism. AqueousKinetics(thermo_t* thermo = 0); - AqueousKinetics(const AqueousKinetics& right); - - AqueousKinetics& operator=(const AqueousKinetics& right); - //! Duplication routine for objects which inherit from Kinetics /*! * This virtual routine can be used to duplicate %Kinetics objects @@ -59,117 +37,25 @@ public: * m_thermo vector within this object */ virtual Kinetics* duplMyselfAsKinetics(const std::vector & tpVector) const; - //@} virtual int type() const { return cAqueousKinetics; } - virtual doublereal reactantStoichCoeff(size_t k, size_t i) const { - return getValue(m_rrxn[k], i, 0.0); - } - - virtual doublereal productStoichCoeff(size_t k, size_t i) const { - return getValue(m_prxn[k], i, 0.0); - } - - //! @name Reaction Rates Of Progress - //@{ virtual void getEquilibriumConstants(doublereal* kc); - - virtual void getDeltaGibbs(doublereal* deltaG); - virtual void getDeltaEnthalpy(doublereal* deltaH); - virtual void getDeltaEntropy(doublereal* deltaS); - - virtual void getDeltaSSGibbs(doublereal* deltaG); - virtual void getDeltaSSEnthalpy(doublereal* deltaH); - virtual void getDeltaSSEntropy(doublereal* deltaS); - - //! @} - //! @name Reaction Mechanism Informational Query Routines - //! @{ - - virtual bool isReversible(size_t i) { - if (std::find(m_revindex.begin(), m_revindex.end(), i) - < m_revindex.end()) { - return true; - } else { - return false; - } - } - virtual void getFwdRateConstants(doublereal* kfwd); - virtual void getRevRateConstants(doublereal* krev, - bool doIrreversible = false); - - //! @} - //! @name Reaction Mechanism Setup Routines - //! @{ - - virtual void init(); - virtual void addReaction(ReactionData& r); - virtual void finalize(); - virtual bool ready() const; - - virtual void update_T(); - - virtual void update_C(); - void updateROP(); - /*! - * Update temperature-dependent portions of reaction rates and - * falloff functions. - */ + //! Update temperature-dependent portions of reaction rates void _update_rates_T(); - /*! - * Update properties that depend on concentrations. Currently only - * the enhanced collision partner concentrations are updated here. - */ + //! Update properties that depend on concentrations. void _update_rates_C(); - //@} - -protected: - - size_t m_nfall; - - Rate1 m_rates; - - std::vector m_irrev; - - size_t m_nirrev; - size_t m_nrev; - - /** - * Difference between the input global reactants order - * and the input global products order. Changed to a double - * to account for the fact that we can have real-valued - * stoichiometries. - */ - vector_fp m_dn; - std::vector m_revindex; - - vector_fp m_conc; - vector_fp m_grt; - - //! @name Aqueous kinetics data - //!@{ - bool m_ROP_ok; - - doublereal m_temp; - //!@} - -private: - void addElementaryReaction(ReactionData& r); - - /** - * Update the equilibrium constants in molar units. - */ + //! Update the equilibrium constants in molar units. void updateKc(); - bool m_finalized; + virtual void addReaction(ReactionData& r); }; } diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h new file mode 100644 index 000000000..935900341 --- /dev/null +++ b/include/cantera/kinetics/BulkKinetics.h @@ -0,0 +1,63 @@ +/** + * @file BulkKinetics.h + * @ingroup chemkinetics + */ + +#ifndef CT_BULKKINETICS_H +#define CT_BULKKINETICS_H + +#include "Kinetics.h" +#include "RateCoeffMgr.h" + +namespace Cantera +{ + +//! Partial specialization of Kinetics for chemistry in a single bulk phase +class BulkKinetics : public Kinetics +{ +public: + BulkKinetics(thermo_t* thermo = 0); + virtual Kinetics* duplMyselfAsKinetics(const std::vector & tpVector) const; + + virtual bool isReversible(size_t i); + + virtual void getDeltaGibbs(doublereal* deltaG); + virtual void getDeltaEnthalpy(doublereal* deltaH); + virtual void getDeltaEntropy(doublereal* deltaS); + + virtual void getDeltaSSGibbs(doublereal* deltaG); + virtual void getDeltaSSEnthalpy(doublereal* deltaH); + virtual void getDeltaSSEntropy(doublereal* deltaS); + + virtual void getRevRateConstants(doublereal* krev, + bool doIrreversible = false); + + virtual void addReaction(ReactionData& r); + virtual void init(); + virtual void finalize(); + virtual bool ready() const; + +protected: + virtual void addElementaryReaction(ReactionData& r); + + Rate1 m_rates; + std::vector m_revindex; //!< Indices of reversible reactions + std::vector m_irrev; //!< Indices of irreversible reactions + + //! Difference between the global reactants order and the global products + //! order. Of type "double" to account for the fact that we can have real- + //! valued stoichiometries. + vector_fp m_dn; + + vector_fp m_conc; + vector_fp m_grt; + + bool m_ROP_ok; + doublereal m_temp; + + bool m_finalized; +}; + +} + +#endif diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index 9ae3ec7ee..4acdcdc41 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -9,26 +9,20 @@ #ifndef CT_GASKINETICS_H #define CT_GASKINETICS_H -#include "Kinetics.h" -#include "ReactionStoichMgr.h" +#include "BulkKinetics.h" #include "ThirdBodyMgr.h" #include "FalloffMgr.h" -#include "RateCoeffMgr.h" namespace Cantera { -// forward references -class Enhanced3BConc; -class ReactionData; - /** * Kinetics manager for elementary gas-phase chemistry. This * kinetics manager implements standard mass-action reaction rate * expressions for low-density gases. * @ingroup kinetics */ -class GasKinetics : public Kinetics +class GasKinetics : public BulkKinetics { public: //! @name Constructors and General Information @@ -40,12 +34,6 @@ public: */ GasKinetics(thermo_t* thermo = 0); - //! Copy Constructor - GasKinetics(const GasKinetics& right); - - //! Assignment operator - GasKinetics& operator=(const GasKinetics& right); - virtual Kinetics* duplMyselfAsKinetics(const std::vector & tpVector) const; virtual int type() const { @@ -57,32 +45,8 @@ public: //! @{ virtual void getEquilibriumConstants(doublereal* kc); - virtual void getDeltaGibbs(doublereal* deltaG); - virtual void getDeltaEnthalpy(doublereal* deltaH); - virtual void getDeltaEntropy(doublereal* deltaS); - - virtual void getDeltaSSGibbs(doublereal* deltaG); - virtual void getDeltaSSEnthalpy(doublereal* deltaH); - virtual void getDeltaSSEntropy(doublereal* deltaS); - - //! @} - //! @name Reaction Mechanism Informational Query Routines - //! @{ - - virtual bool isReversible(size_t i) { - if (std::find(m_revindex.begin(), m_revindex.end(), i) - < m_revindex.end()) { - return true; - } else { - return false; - } - } - virtual void getFwdRateConstants(doublereal* kfwd); - virtual void getRevRateConstants(doublereal* krev, - bool doIrreversible = false); - //! @} //! @name Reaction Mechanism Setup Routines //! @{ @@ -111,30 +75,15 @@ protected: Rate1 m_falloff_low_rates; Rate1 m_falloff_high_rates; - Rate1 m_rates; FalloffMgr m_falloffn; ThirdBodyMgr m_3b_concm; ThirdBodyMgr m_falloff_concm; - std::vector m_irrev; - Rate1 m_plog_rates; Rate1 m_cheb_rates; - size_t m_nirrev; - size_t m_nrev; - - /** - * Difference between the input global reactants order - * and the input global products order. Changed to a double - * to account for the fact that we can have real-valued - * stoichiometries. - */ - vector_fp m_dn; - std::vector m_revindex; - //! @name Reaction rate data //!@{ doublereal m_logp_ref; @@ -142,21 +91,15 @@ protected: doublereal m_logStandConc; vector_fp m_rfn_low; vector_fp m_rfn_high; - bool m_ROP_ok; - doublereal m_temp; doublereal m_pres; //!< Last pressure at which rates were evaluated vector_fp falloff_work; vector_fp concm_3b_values; vector_fp concm_falloff_values; //!@} - vector_fp m_conc; void processFalloffReactions(); - vector_fp m_grt; -private: - void addElementaryReaction(ReactionData& r); void addThreeBodyReaction(ReactionData& r); void addFalloffReaction(ReactionData& r); void addPlogReaction(ReactionData& r); diff --git a/src/kinetics/AqueousKinetics.cpp b/src/kinetics/AqueousKinetics.cpp index 4d86e6ce5..12e0c9669 100644 --- a/src/kinetics/AqueousKinetics.cpp +++ b/src/kinetics/AqueousKinetics.cpp @@ -3,7 +3,6 @@ * * Homogeneous kinetics in an aqueous phase, either condensed * or dilute in salts - * */ /* * Copyright (2006) Sandia Corporation. Under the terms of @@ -20,58 +19,8 @@ namespace Cantera { AqueousKinetics::AqueousKinetics(thermo_t* thermo) : - m_nfall(0), - m_nirrev(0), - m_nrev(0), - m_ROP_ok(false), - m_temp(0.0), - m_finalized(false) + BulkKinetics(thermo) { - if (thermo != 0) { - addPhase(*thermo); - } -} - -AqueousKinetics::AqueousKinetics(const AqueousKinetics& right) : - m_nfall(0), - m_nirrev(0), - m_nrev(0), - m_ROP_ok(false), - m_temp(0.0), - m_finalized(false) -{ - *this = right; -} - -AqueousKinetics& AqueousKinetics::operator=(const AqueousKinetics& right) -{ - if (this == &right) { - return *this; - } - - Kinetics::operator=(right); - - m_nfall = right.m_nfall; - m_rates = right.m_rates; - m_irrev = right.m_irrev; - m_nirrev = right.m_nirrev; - m_nrev = right.m_nrev; - m_rxntype = right.m_rxntype; - m_dn = right.m_dn; - m_revindex = right.m_revindex; - - m_ROP_ok = right.m_ROP_ok; - m_temp = right.m_temp; - - m_conc = right.m_conc; - m_grt = right.m_grt; - m_finalized = right.m_finalized; - - throw CanteraError("GasKinetics::operator=()", - "Unfinished implementation"); - - return *this; - } Kinetics* AqueousKinetics::duplMyselfAsKinetics(const std::vector & tpVector) const @@ -81,19 +30,10 @@ Kinetics* AqueousKinetics::duplMyselfAsKinetics(const std::vector & t return gK; } -void AqueousKinetics::update_T() -{ -} - -void AqueousKinetics::update_C() -{ -} - void AqueousKinetics::_update_rates_T() { doublereal T = thermo().temperature(); - doublereal logT = log(T); - m_rates.update(T, logT, &m_rfn[0]); + m_rates.update(T, log(T), &m_rfn[0]); m_temp = T; updateKc(); @@ -103,7 +43,6 @@ void AqueousKinetics::_update_rates_T() void AqueousKinetics::_update_rates_C() { thermo().getActivityConcentrations(&m_conc[0]); - m_ROP_ok = false; } @@ -122,12 +61,12 @@ void AqueousKinetics::updateKc() m_rxnstoich.getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]); doublereal rrt = 1.0/(GasConstant * thermo().temperature()); - for (size_t i = 0; i < m_nrev; i++) { + for (size_t i = 0; i < m_revindex.size(); i++) { size_t irxn = m_revindex[i]; m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt); } - for (size_t i = 0; i != m_nirrev; ++i) { + for (size_t i = 0; i != m_irrev.size(); ++i) { m_rkcn[ m_irrev[i] ] = 0.0; } } @@ -157,103 +96,6 @@ void AqueousKinetics::getEquilibriumConstants(doublereal* kc) m_temp = 0.0; } -void AqueousKinetics::getDeltaGibbs(doublereal* deltaG) -{ - /* - * Get the chemical potentials of the species in the - * ideal gas solution. - */ - thermo().getChemPotentials(&m_grt[0]); - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG); -} - -void AqueousKinetics::getDeltaEnthalpy(doublereal* deltaH) -{ - /* - * Get the partial molar enthalpy of all species in the - * ideal gas. - */ - thermo().getPartialMolarEnthalpies(&m_grt[0]); - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH); -} - -void AqueousKinetics::getDeltaEntropy(doublereal* deltaS) -{ - /* - * Get the partial molar entropy of all species in the - * solid solution. - */ - thermo().getPartialMolarEntropies(&m_grt[0]); - /* - * Use the stoichiometric manager to find deltaS for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS); -} - -void AqueousKinetics::getDeltaSSGibbs(doublereal* deltaG) -{ - /* - * Get the standard state chemical potentials of the species. - * This is the array of chemical potentials at unit activity - * We define these here as the chemical potentials of the pure - * species at the temperature and pressure of the solution. - */ - thermo().getStandardChemPotentials(&m_grt[0]); - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG); -} - -void AqueousKinetics::getDeltaSSEnthalpy(doublereal* deltaH) -{ - /* - * Get the standard state enthalpies of the species. - * This is the array of chemical potentials at unit activity - * We define these here as the enthalpies of the pure - * species at the temperature and pressure of the solution. - */ - thermo().getEnthalpy_RT(&m_grt[0]); - doublereal RT = thermo().temperature() * GasConstant; - for (size_t k = 0; k < m_kk; k++) { - m_grt[k] *= RT; - } - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH); -} - -void AqueousKinetics::getDeltaSSEntropy(doublereal* deltaS) -{ - /* - * Get the standard state entropy of the species. - * We define these here as the entropies of the pure - * species at the temperature and pressure of the solution. - */ - thermo().getEntropy_R(&m_grt[0]); - doublereal R = GasConstant; - for (size_t k = 0; k < m_kk; k++) { - m_grt[k] *= R; - } - /* - * Use the stoichiometric manager to find deltaS for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS); -} - void AqueousKinetics::updateROP() { _update_rates_T(); @@ -272,19 +114,15 @@ void AqueousKinetics::updateROP() // copy the forward rates to the reverse rates copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin()); - // for reverse rates computed from thermochemistry, multiply - // the forward rates copied into m_ropr by the reciprocals of - // the equilibrium constants + // for reverse rates computed from thermochemistry, multiply the forward + // rates copied into m_ropr by the reciprocals of the equilibrium constants multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin()); // multiply ropf by concentration products m_rxnstoich.multiplyReactants(&m_conc[0], &m_ropf[0]); - //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); - // for reversible reactions, multiply ropr by concentration - // products + // for reversible reactions, multiply ropr by concentration products m_rxnstoich.multiplyRevProducts(&m_conc[0], &m_ropr[0]); - //m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); for (size_t j = 0; j != m_ii; ++j) { m_ropnet[j] = m_ropf[j] - m_ropr[j]; @@ -309,84 +147,13 @@ void AqueousKinetics::getFwdRateConstants(doublereal* kfwd) } } -void AqueousKinetics::getRevRateConstants(doublereal* krev, - bool doIrreversible) -{ - /* - * go get the forward rate constants. -> note, we don't - * really care about speed or redundancy in these - * informational routines. - */ - getFwdRateConstants(krev); - - if (doIrreversible) { - getEquilibriumConstants(&m_ropnet[0]); - for (size_t i = 0; i < m_ii; i++) { - krev[i] /= m_ropnet[i]; - } - } else { - /* - * m_rkcn[] is zero for irreversible reactions - */ - for (size_t i = 0; i < m_ii; i++) { - krev[i] *= m_rkcn[i]; - } - } -} - void AqueousKinetics::addReaction(ReactionData& r) { if (r.reactionType == ELEMENTARY_RXN) { addElementaryReaction(r); } - m_dn.push_back(accumulate(r.pstoich.begin(), r.pstoich.end(), 0.0) - - accumulate(r.rstoich.begin(), r.rstoich.end(), 0.0)); - if (r.reversible) { - m_revindex.push_back(nReactions()); - m_nrev++; - } else { - m_irrev.push_back(nReactions()); - m_nirrev++; - } - Kinetics::addReaction(r); -} - -void AqueousKinetics::addElementaryReaction(ReactionData& r) -{ - // install rate coeff calculator - m_rates.install(nReactions(), r); -} - -void AqueousKinetics::init() -{ - m_kk = thermo().nSpecies(); - m_rrxn.resize(m_kk); - m_prxn.resize(m_kk); - m_conc.resize(m_kk); - m_grt.resize(m_kk); -} - -void AqueousKinetics::finalize() -{ - if (!m_finalized) { - m_finalized = true; - - // Guarantee that these arrays can be converted to double* even in the - // special case where there are no reactions defined. - if (!m_ii) { - m_perturb.resize(1, 1.0); - m_ropf.resize(1, 0.0); - m_ropr.resize(1, 0.0); - m_ropnet.resize(1, 0.0); - m_rkcn.resize(1, 0.0); - } - } -} - -bool AqueousKinetics::ready() const -{ - return m_finalized; + BulkKinetics::addReaction(r); } } diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp new file mode 100644 index 000000000..edee5c7dd --- /dev/null +++ b/src/kinetics/BulkKinetics.cpp @@ -0,0 +1,164 @@ +#include "cantera/kinetics/BulkKinetics.h" + +namespace Cantera +{ + +BulkKinetics::BulkKinetics(thermo_t* thermo) : + m_ROP_ok(false), + m_temp(0.0), + m_finalized(false) +{ + if (thermo) { + addPhase(*thermo); + } +} + + +Kinetics* BulkKinetics::duplMyselfAsKinetics(const std::vector & tpVector) const +{ + BulkKinetics* kin = new BulkKinetics(*this); + kin->assignShallowPointers(tpVector); + return kin; +} + +bool BulkKinetics::isReversible(size_t i) { + if (std::find(m_revindex.begin(), m_revindex.end(), i) + < m_revindex.end()) { + return true; + } else { + return false; + } +} + + +void BulkKinetics::getDeltaGibbs(doublereal* deltaG) +{ + // Get the chemical potentials of the species in the ideal gas solution. + thermo().getChemPotentials(&m_grt[0]); + // Use the stoichiometric manager to find deltaG for each reaction. + m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG); +} + +void BulkKinetics::getDeltaEnthalpy(doublereal* deltaH) +{ + // Get the partial molar enthalpy of all species in the ideal gas. + thermo().getPartialMolarEnthalpies(&m_grt[0]); + // Use the stoichiometric manager to find deltaH for each reaction. + m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH); +} + +void BulkKinetics::getDeltaEntropy(doublereal* deltaS) +{ + // Get the partial molar entropy of all species in the solid solution. + thermo().getPartialMolarEntropies(&m_grt[0]); + // Use the stoichiometric manager to find deltaS for each reaction. + m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS); +} + +void BulkKinetics::getDeltaSSGibbs(doublereal* deltaG) +{ + // Get the standard state chemical potentials of the species. This is the + // array of chemical potentials at unit activity. We define these here as + // the chemical potentials of the pure species at the temperature and + // pressure of the solution. + thermo().getStandardChemPotentials(&m_grt[0]); + // Use the stoichiometric manager to find deltaG for each reaction. + m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG); +} + +void BulkKinetics::getDeltaSSEnthalpy(doublereal* deltaH) +{ + // Get the standard state enthalpies of the species. + thermo().getEnthalpy_RT(&m_grt[0]); + doublereal RT = thermo().temperature() * GasConstant; + for (size_t k = 0; k < m_kk; k++) { + m_grt[k] *= RT; + } + // Use the stoichiometric manager to find deltaH for each reaction. + m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH); +} + +void BulkKinetics::getDeltaSSEntropy(doublereal* deltaS) +{ + // Get the standard state entropy of the species. We define these here as + // the entropies of the pure species at the temperature and pressure of the + // solution. + thermo().getEntropy_R(&m_grt[0]); + doublereal R = GasConstant; + for (size_t k = 0; k < m_kk; k++) { + m_grt[k] *= R; + } + // Use the stoichiometric manager to find deltaS for each reaction. + m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS); +} + +void BulkKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible) +{ + /* + * go get the forward rate constants. -> note, we don't + * really care about speed or redundancy in these + * informational routines. + */ + getFwdRateConstants(krev); + + if (doIrreversible) { + getEquilibriumConstants(&m_ropnet[0]); + for (size_t i = 0; i < m_ii; i++) { + krev[i] /= m_ropnet[i]; + } + } else { + // m_rkcn[] is zero for irreversible reactions + for (size_t i = 0; i < m_ii; i++) { + krev[i] *= m_rkcn[i]; + } + } +} + +void BulkKinetics::addReaction(ReactionData& r) +{ + m_dn.push_back(accumulate(r.pstoich.begin(), r.pstoich.end(), 0.0) - + accumulate(r.rstoich.begin(), r.rstoich.end(), 0.0)); + + if (r.reversible) { + m_revindex.push_back(nReactions()); + } else { + m_irrev.push_back(nReactions()); + } + Kinetics::addReaction(r); +} + +void BulkKinetics::addElementaryReaction(ReactionData& r) +{ + m_rates.install(nReactions(), r); +} + +void BulkKinetics::init() +{ + m_kk = thermo().nSpecies(); + m_rrxn.resize(m_kk); + m_prxn.resize(m_kk); + m_conc.resize(m_kk); + m_grt.resize(m_kk); +} + +void BulkKinetics::finalize() +{ + m_finalized = true; + + // Guarantee that these arrays can be converted to double* even in the + // special case where there are no reactions defined. + if (!m_ii) { + m_perturb.resize(1, 1.0); + m_ropf.resize(1, 0.0); + m_ropr.resize(1, 0.0); + m_ropnet.resize(1, 0.0); + m_rkcn.resize(1, 0.0); + } +} + +bool BulkKinetics::ready() const +{ + return m_finalized; +} + +} diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index a379161cd..21222e3e6 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -13,73 +13,13 @@ using namespace std; namespace Cantera { GasKinetics::GasKinetics(thermo_t* thermo) : + BulkKinetics(thermo), m_nfall(0), - m_nirrev(0), - m_nrev(0), m_logp_ref(0.0), m_logc_ref(0.0), m_logStandConc(0.0), - m_ROP_ok(false), - m_temp(0.0), - m_pres(0.0), - m_finalized(false) + m_pres(0.0) { - if (thermo != 0) { - addPhase(*thermo); - } -} - -GasKinetics::GasKinetics(const GasKinetics& right) -{ - *this = right; -} - -GasKinetics& GasKinetics::operator=(const GasKinetics& right) -{ - if (this == &right) { - return *this; - } - - Kinetics::operator=(right); - - m_nfall = right.m_nfall; - m_fallindx = right.m_fallindx; - m_falloff_low_rates = right.m_falloff_low_rates; - m_falloff_high_rates = right.m_falloff_high_rates; - m_rates = right.m_rates; - m_falloffn = right.m_falloffn; - m_3b_concm = right.m_3b_concm; - m_falloff_concm = right.m_falloff_concm; - m_irrev = right.m_irrev; - m_plog_rates = right.m_plog_rates; - m_cheb_rates = right.m_cheb_rates; - - m_nirrev = right.m_nirrev; - m_nrev = right.m_nrev; - m_rrxn = right.m_rrxn; - m_prxn = right.m_prxn; - m_dn = right.m_dn; - m_revindex = right.m_revindex; - - m_logp_ref = right.m_logp_ref; - m_logc_ref = right.m_logc_ref; - m_logStandConc = right.m_logStandConc; - m_rfn_low = right.m_rfn_low; - m_rfn_high = right.m_rfn_high; - m_ROP_ok = right.m_ROP_ok; - m_temp = right.m_temp; - falloff_work = right.falloff_work; - concm_3b_values = right.concm_3b_values; - concm_falloff_values = right.concm_falloff_values; - - m_conc = right.m_conc; - m_grt = right.m_grt; - m_finalized = right.m_finalized; - - throw CanteraError("GasKinetics::operator=()", - "Unfinished implementation"); - - return *this; } Kinetics* GasKinetics::duplMyselfAsKinetics(const std::vector & tpVector) const @@ -166,13 +106,13 @@ void GasKinetics::updateKc() m_rxnstoich.getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]); doublereal rrt = 1.0/(GasConstant * thermo().temperature()); - for (size_t i = 0; i < m_nrev; i++) { + for (size_t i = 0; i < m_revindex.size(); i++) { size_t irxn = m_revindex[i]; m_rkcn[irxn] = std::min(exp(m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc), BigNumber); } - for (size_t i = 0; i != m_nirrev; ++i) { + for (size_t i = 0; i != m_irrev.size(); ++i) { m_rkcn[ m_irrev[i] ] = 0.0; } } @@ -196,103 +136,6 @@ void GasKinetics::getEquilibriumConstants(doublereal* kc) m_temp = 0.0; } -void GasKinetics::getDeltaGibbs(doublereal* deltaG) -{ - /* - * Get the chemical potentials of the species in the - * ideal gas solution. - */ - thermo().getChemPotentials(&m_grt[0]); - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG); -} - -void GasKinetics::getDeltaEnthalpy(doublereal* deltaH) -{ - /* - * Get the partial molar enthalpy of all species in the - * ideal gas. - */ - thermo().getPartialMolarEnthalpies(&m_grt[0]); - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH); -} - -void GasKinetics::getDeltaEntropy(doublereal* deltaS) -{ - /* - * Get the partial molar entropy of all species in the - * solid solution. - */ - thermo().getPartialMolarEntropies(&m_grt[0]); - /* - * Use the stoichiometric manager to find deltaS for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS); -} - -void GasKinetics::getDeltaSSGibbs(doublereal* deltaG) -{ - /* - * Get the standard state chemical potentials of the species. - * This is the array of chemical potentials at unit activity - * We define these here as the chemical potentials of the pure - * species at the temperature and pressure of the solution. - */ - thermo().getStandardChemPotentials(&m_grt[0]); - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG); -} - -void GasKinetics::getDeltaSSEnthalpy(doublereal* deltaH) -{ - /* - * Get the standard state enthalpies of the species. - * This is the array of chemical potentials at unit activity - * We define these here as the enthalpies of the pure - * species at the temperature and pressure of the solution. - */ - thermo().getEnthalpy_RT(&m_grt[0]); - doublereal RT = thermo().temperature() * GasConstant; - for (size_t k = 0; k < m_kk; k++) { - m_grt[k] *= RT; - } - /* - * Use the stoichiometric manager to find deltaG for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH); -} - -void GasKinetics::getDeltaSSEntropy(doublereal* deltaS) -{ - /* - * Get the standard state entropy of the species. - * We define these here as the entropies of the pure - * species at the temperature and pressure of the solution. - */ - thermo().getEntropy_R(&m_grt[0]); - doublereal R = GasConstant; - for (size_t k = 0; k < m_kk; k++) { - m_grt[k] *= R; - } - /* - * Use the stoichiometric manager to find deltaS for each - * reaction. - */ - m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS); -} - void GasKinetics::processFalloffReactions() { // use m_ropr for temporary storage of reduced pressure @@ -346,19 +189,15 @@ void GasKinetics::updateROP() // copy the forward rates to the reverse rates copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin()); - // for reverse rates computed from thermochemistry, multiply - // the forward rates copied into m_ropr by the reciprocals of - // the equilibrium constants + // for reverse rates computed from thermochemistry, multiply the forward + // rates copied into m_ropr by the reciprocals of the equilibrium constants multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin()); // multiply ropf by concentration products m_rxnstoich.multiplyReactants(&m_conc[0], &m_ropf[0]); - //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); - // for reversible reactions, multiply ropr by concentration - // products + // for reversible reactions, multiply ropr by concentration products m_rxnstoich.multiplyRevProducts(&m_conc[0], &m_ropr[0]); - //m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); for (size_t j = 0; j != m_ii; ++j) { m_ropnet[j] = m_ropf[j] - m_ropr[j]; @@ -389,10 +228,6 @@ void GasKinetics::getFwdRateConstants(doublereal* kfwd) m_3b_concm.multiply(&m_ropf[0], &concm_3b_values[0]); } - /* - * This routine is hardcoded to replace some of the values - * of the ropf vector. - */ if (m_nfall) { processFalloffReactions(); } @@ -405,28 +240,6 @@ void GasKinetics::getFwdRateConstants(doublereal* kfwd) } } -void GasKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible) -{ - /* - * go get the forward rate constants. -> note, we don't - * really care about speed or redundancy in these - * informational routines. - */ - getFwdRateConstants(krev); - - if (doIrreversible) { - getEquilibriumConstants(&m_ropnet[0]); - for (size_t i = 0; i < m_ii; i++) { - krev[i] /= m_ropnet[i]; - } - } else { - // m_rkcn[] is zero for irreversible reactions - for (size_t i = 0; i < m_ii; i++) { - krev[i] *= m_rkcn[i]; - } - } -} - void GasKinetics::addReaction(ReactionData& r) { switch (r.reactionType) { @@ -451,17 +264,7 @@ void GasKinetics::addReaction(ReactionData& r) } // operations common to all reaction types - m_dn.push_back(accumulate(r.pstoich.begin(), r.pstoich.end(), 0.0) - - accumulate(r.rstoich.begin(), r.rstoich.end(), 0.0)); - - if (r.reversible) { - m_revindex.push_back(nReactions()); - m_nrev++; - } else { - m_irrev.push_back(nReactions()); - m_nirrev++; - } - Kinetics::addReaction(r); + BulkKinetics::addReaction(r); } void GasKinetics::addFalloffReaction(ReactionData& r) @@ -474,17 +277,15 @@ void GasKinetics::addFalloffReaction(ReactionData& r) m_falloff_low_rates.install(m_nfall, r); m_rfn_low.push_back(r.rateCoeffParameters[0]); - // add this reaction number to the list of - // falloff reactions + // add this reaction number to the list of falloff reactions m_fallindx.push_back(nReactions()); - // install the enhanced third-body concentration - // calculator for this reaction + // install the enhanced third-body concentration calculator for this + // reaction m_falloff_concm.install(m_nfall, r.thirdBodyEfficiencies, r.default_3b_eff); - // install the falloff function calculator for - // this reaction + // install the falloff function calculator for this reaction m_falloffn.install(m_nfall, r.falloffType, r.reactionType, r.falloffParameters); @@ -492,14 +293,8 @@ void GasKinetics::addFalloffReaction(ReactionData& r) ++m_nfall; } -void GasKinetics::addElementaryReaction(ReactionData& r) -{ - m_rates.install(nReactions(), r); -} - void GasKinetics::addThreeBodyReaction(ReactionData& r) { - // install rate coeff calculator m_rates.install(nReactions(), r); m_3b_concm.install(nReactions(), r.thirdBodyEfficiencies, r.default_3b_eff); @@ -507,44 +302,26 @@ void GasKinetics::addThreeBodyReaction(ReactionData& r) void GasKinetics::addPlogReaction(ReactionData& r) { - // install rate coefficient calculator m_plog_rates.install(nReactions(), r); } void GasKinetics::addChebyshevReaction(ReactionData& r) { - // install rate coefficient calculator m_cheb_rates.install(nReactions(), r); } void GasKinetics::init() { - m_kk = thermo().nSpecies(); - m_rrxn.resize(m_kk); - m_prxn.resize(m_kk); - m_conc.resize(m_kk); - m_grt.resize(m_kk); + BulkKinetics::init(); m_logp_ref = log(thermo().refPressure()) - log(GasConstant); } void GasKinetics::finalize() { - if (!m_finalized) { - falloff_work.resize(m_falloffn.workSize()); - concm_3b_values.resize(m_3b_concm.workSize()); - concm_falloff_values.resize(m_falloff_concm.workSize()); - m_finalized = true; - - // Guarantee that these arrays can be converted to double* even in the - // special case where there are no reactions defined. - if (!m_ii) { - m_perturb.resize(1, 1.0); - m_ropf.resize(1, 0.0); - m_ropr.resize(1, 0.0); - m_ropnet.resize(1, 0.0); - m_rkcn.resize(1, 0.0); - } - } + BulkKinetics::finalize(); + falloff_work.resize(m_falloffn.workSize()); + concm_3b_values.resize(m_3b_concm.workSize()); + concm_falloff_values.resize(m_falloff_concm.workSize()); } bool GasKinetics::ready() const