From 510b25b884e3366200310979eac2dc0df8637ac9 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 30 Mar 2012 23:47:50 +0000 Subject: [PATCH] Implementation of P-log rate expressions Also includes skeleton for implementing Chebyshev rates. --- include/cantera/kinetics/GasKinetics.h | 8 + include/cantera/kinetics/RateCoeffMgr.h | 4 + include/cantera/kinetics/ReactionData.h | 4 +- include/cantera/kinetics/RxnRates.h | 233 ++++++++++++++++++++++++ src/kinetics/GasKinetics.cpp | 58 ++++-- src/kinetics/importKinetics.cpp | 3 +- 6 files changed, 297 insertions(+), 13 deletions(-) diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index 90bbcee05..dc85d33ad 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -386,6 +386,11 @@ public: void _update_rates_T(); + + //! Update properties that depend on concentrations. + //! Currently the enhanced collision partner concentrations are updated + //! here, as well as the pressure-dependent portion of P-log and Chebyshev + //! reactions. void _update_rates_C(); //@} @@ -409,6 +414,9 @@ protected: std::vector m_irrev; + Rate1 m_plog_rates; + Rate1 m_cheb_rates; + ReactionStoichMgr* m_rxnstoich; std::vector m_fwdOrder; diff --git a/include/cantera/kinetics/RateCoeffMgr.h b/include/cantera/kinetics/RateCoeffMgr.h index 6646bc984..6bc1df8cc 100644 --- a/include/cantera/kinetics/RateCoeffMgr.h +++ b/include/cantera/kinetics/RateCoeffMgr.h @@ -104,6 +104,10 @@ public: output1 << key; } + size_t nReactions() const { + return m_rates.size(); + } + protected: std::vector m_rates; std::vector m_rxn; diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h index 770e10d8e..cc58693c6 100644 --- a/include/cantera/kinetics/ReactionData.h +++ b/include/cantera/kinetics/ReactionData.h @@ -74,7 +74,9 @@ public: //! Arrhenius parameters for P-log reactions. //! The keys are the pressures corresponding to each Arrhenius expression. - std::map plogParameters; + //! Multiple sets of Arrhenius parameters may be specified at a given + //! pressure. + std::multimap plogParameters; double chebTmin; //!< Minimum temperature for Chebyshev fit double chebTmax; //!< Maximum temperature for Chebyshev fit diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index bb401fbb5..c1872b748 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -468,6 +468,239 @@ protected: }; +class Plog +{ +public: + //! return the rate coefficient type. + static int type() + { + return PLOG_REACTION_RATECOEFF_TYPE; + } + + //! Default constructor. + Plog() {} + + //! Constructor from ReactionData. + explicit Plog(const ReactionData& rdata) : + logP1_(1000), + logP2_(-1000), + maxRates_(1) + { + typedef std::multimap::const_iterator iter_t; + + size_t j = 0; + size_t rateCount = 0; + // Insert intermediate pressures + for (iter_t iter = rdata.plogParameters.begin(); + iter != rdata.plogParameters.end(); + iter++) { + double logp = log(iter->first); + if (pressures_.empty() || pressures_.rbegin()->first != logp) { + // starting a new group + pressures_[logp] = std::make_pair(j, j+1); + rateCount = 1; + } else { + // another rate expression at the same pressure + pressures_[logp].second = j+1; + rateCount++; + } + maxRates_ = std::max(rateCount, maxRates_); + + j++; + A_.push_back(iter->second[0]); + n_.push_back(iter->second[1]); + Ea_.push_back(iter->second[2]); + } + + // For pressures with only one Arrhenius expression, it is more + // efficient to work with log(A) + for (pressureIter iter = pressures_.begin(); + iter != pressures_.end(); + iter++) { + if (iter->second.first == iter->second.second - 1) { + A_[iter->second.first] = log(A_[iter->second.first]); + } + } + + // Duplicate the first and last groups to handle P < P_0 and P > P_N + pressures_.insert(std::make_pair(-1000.0, pressures_.begin()->second)); + pressures_.insert(std::make_pair(1000.0, pressures_.rbegin()->second)); + + // Resize work arrays + A1_.resize(maxRates_); + A2_.resize(maxRates_); + n1_.resize(maxRates_); + n2_.resize(maxRates_); + Ea1_.resize(maxRates_); + Ea2_.resize(maxRates_); + } + + //! 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; + } + + pressureIter iter = pressures_.upper_bound(c[0]); + AssertThrowMsg(iter != pressures_.end(), "Plog::update_C", + "Pressure out of range: " + fp2str(logP)); + AssertThrowMsg(iter != pressures.begin(), "Plog::update_C", + "Pressure out of range: " + fp2str(logP)); + + // upper interpolation pressure + logP2_ = iter->first; + size_t start = iter->second.first; + m2_ = iter->second.second - start; + for (size_t m = 0; m < m2_; m++) { + A2_[m] = A_[start+m]; + n2_[m] = n_[start+m]; + Ea2_[m] = Ea_[start+m]; + } + + // lower interpolation pressure + logP1_ = (--iter)->first; + start = iter->second.first; + m1_ = iter->second.second - start; + for (size_t m = 0; m < m1_; m++) { + A1_[m] = A_[start+m]; + n1_[m] = n_[start+m]; + Ea1_[m] = Ea_[start+m]; + } + + rDeltaP_ = 1.0 / (logP2_ - logP1_); + } + + /** + * Update the value of the logarithm of the rate constant. + */ + doublereal update(doublereal logT, doublereal recipT) const + { + double log_k1, log_k2; + if (m1_ == 1) { + log_k1 = A1_[0] + n1_[0] * logT - Ea1_[0] * recipT; + } else { + double k = 0.0; + for (size_t m = 0; m < m1_; m++) { + k += A1_[m] * exp(n1_[m] * logT - Ea1_[m] * recipT); + } + log_k1 = log(k); + } + + if (m2_ == 1) { + log_k2 = A2_[0] + n2_[0] * logT - Ea2_[0] * recipT; + } else { + double k = 0.0; + for (size_t m = 0; m < m2_; m++) { + k += A2_[m] * exp(n2_[m] * logT - Ea2_[m] * recipT); + } + log_k2 = log(k); + } + + return log_k1 + (log_k2 - log_k1) * (logP_ - logP1_) * rDeltaP_; + } + + /** + * Update the value the rate constant. + * + * This function returns the actual value of the rate constant. + */ + doublereal updateRC(doublereal logT, doublereal recipT) const { + return exp(update(logT, recipT)); + } + + doublereal activationEnergy_R() const { + throw CanteraError("Plog::activationEnergy_R", "Not implemented"); + } + + static bool alwaysComputeRate() { + return false; + } + +protected: + //! log(p) to (index range) in A_, n, Ea vectors + std::map > pressures_; + typedef std::map >::iterator pressureIter; + + vector_fp A_; //!< Pre-exponential factor at each pressure (or log(A)) + vector_fp n_; //!< Temperature exponent at each pressure [dimensionless] + vector_fp Ea_; //!< Activation energy at each pressure [K] + + double logP_; //!< log(p) at the current state + double logP1_, logP2_; //!< log(p) at the lower / upper pressure reference + + //! Pre-exponential factors at lower / upper pressure reference. + //! Stored as log(A) when there is only one at the corresponding pressure. + vector_fp A1_, A2_; + vector_fp n1_, n2_; //!< n at lower / upper pressure reference + vector_fp Ea1_, Ea2_; //!< Activation energy at lower / upper pressure reference + + //! Number of Arrhenius expressions at lower / upper pressure references + size_t m1_, m2_; + double rDeltaP_; //!< reciprocal of (logP2 - logP1) + + size_t maxRates_; //!< The maximum number of rates at any given pressure +}; + + +class ChebyshevRate +{ +public: + //! return the rate coefficient type. + static int type() + { + return CHEBYSHEV_REACTION_RATECOEFF_TYPE; + } + + //! Default constructor. + ChebyshevRate() {} + + //! Constructor from ReactionData. + explicit ChebyshevRate(const ReactionData& rdata) + { + } + + //! Update concentration-dependent parts of the rate coefficient. + //! @param c natural log of the pressure in Pa + void update_C(const doublereal* c) + { + } + + /** + * Update the value of the logarithm of the rate constant. + * + * Note, this function should never be called for negative A values. + * If it does then it will produce a negative overflow result, and + * a zero net forwards reaction rate, instead of a negative reaction + * rate constant that is the expected result. + */ + doublereal update(doublereal logT, doublereal recipT) const + { + return 0.0; + } + + /** + * Update the value the rate constant. + * + * This function returns the actual value of the rate constant. + */ + doublereal updateRC(doublereal logT, doublereal recipT) const { + return exp(update(logT, recipT)); + } + + doublereal activationEnergy_R() const { + return 0.0; + } + + static bool alwaysComputeRate() { + return false; + } + +protected: +}; + // class LandauTeller { // public: diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 590acf1ec..52acf39f5 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -129,6 +129,8 @@ GasKinetics& GasKinetics::operator=(const GasKinetics& right) 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_rxnstoich = *(right.m_rxnstoich); @@ -195,6 +197,7 @@ _update_rates_T() if (!m_kdata->m_rfn.empty()) { m_rates.update(T, logT, &m_kdata->m_rfn[0]); } + if (!m_kdata->m_rfn_low.empty()) { m_falloff_low_rates.update(T, logT, &m_kdata->m_rfn_low[0]); m_falloff_high_rates.update(T, logT, &m_kdata->m_rfn_high[0]); @@ -202,29 +205,50 @@ _update_rates_T() if (!m_kdata->falloff_work.empty()) { m_falloffn.updateTemp(T, &m_kdata->falloff_work[0]); } + if (m_plog_rates.nReactions()) { + m_plog_rates.update(T, logT, &m_kdata->m_rfn[0]); + } + + if (m_cheb_rates.nReactions()) { + m_cheb_rates.update(T, logT, &m_kdata->m_rfn[0]); + } + m_kdata->m_temp = T; updateKc(); m_kdata->m_ROP_ok = false; - //} }; //==================================================================================================================== -/** - * Update properties that depend on concentrations. Currently only - * the enhanced collision partner concentrations are updated here. - */ + void GasKinetics:: _update_rates_C() { thermo().getActivityConcentrations(&m_conc[0]); doublereal ctot = thermo().molarDensity(); + + // 3-body reactions if (!m_kdata->concm_3b_values.empty()) { m_3b_concm.update(m_conc, ctot, &m_kdata->concm_3b_values[0]); } + + // Falloff reactions if (!m_kdata->concm_falloff_values.empty()) { m_falloff_concm.update(m_conc, ctot, &m_kdata->concm_falloff_values[0]); } + + double logP = log(thermo().pressure()); + + // P-log reactions + if (m_plog_rates.nReactions()) { + m_plog_rates.update_C(&logP); + } + + // Chebyshev reactions + if (m_cheb_rates.nReactions()) { + m_cheb_rates.update_C(&logP); + } + m_kdata->m_ROP_ok = false; } //==================================================================================================================== @@ -523,9 +547,8 @@ void GasKinetics::processFalloffReactions() //==================================================================================================================== void GasKinetics::updateROP() { - - _update_rates_T(); _update_rates_C(); + _update_rates_T(); if (m_kdata->m_ROP_ok) { return; @@ -587,8 +610,8 @@ void GasKinetics::updateROP() void GasKinetics:: getFwdRateConstants(doublereal* kfwd) { - _update_rates_T(); _update_rates_C(); + _update_rates_T(); // copy rate coefficients into ropf const vector_fp& rf = m_kdata->m_rfn; @@ -761,17 +784,30 @@ addThreeBodyReaction(ReactionData& r) void GasKinetics::addPlogReaction(ReactionData& r) { - // @todo: Not yet implemented + // install rate coefficient calculator + size_t iloc = m_plog_rates.install(reactionNumber(), r); + + // add a dummy entry in m_rfn, where computed rate coeff will be put + m_kdata->m_rfn.push_back(0.0); + + m_fwdOrder.push_back(r.reactants.size()); + registerReaction(reactionNumber(), PLOG_RXN, iloc); } void GasKinetics::addChebyshevReaction(ReactionData& r) { - // @todo: Not yet implemented + // install rate coefficient calculator + size_t iloc = m_cheb_rates.install(reactionNumber(), r); + + // add a dummy entry in m_rfn, where computed rate coeff will be put + m_kdata->m_rfn.push_back(0.0); + + m_fwdOrder.push_back(r.reactants.size()); + registerReaction(reactionNumber(), CHEBYSHEV_RXN, iloc); } void GasKinetics::installReagents(const ReactionData& r) { - m_kdata->m_ropf.push_back(0.0); // extend by one for new rxn m_kdata->m_ropr.push_back(0.0); m_kdata->m_ropnet.push_back(0.0); diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index 25d48fd66..563ed1d4c 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -509,7 +509,8 @@ void getRateCoefficient(const XML_Node& kf, Kinetics& kin, for (size_t m = 0; m < kf.nChildren(); m++) { const XML_Node& node = kf.child(m); double p = getFloat(node, "P", "toSI"); - vector_fp& rate = rdata.plogParameters[p]; + vector_fp& rate = rdata.plogParameters.insert( + std::make_pair(p, vector_fp()))->second; rate.resize(3); rate[0] = getFloat(node, "A", "toSI"); rate[1] = getFloat(node, "b");