diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 8044e6d28..e49ad18d9 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -71,6 +71,13 @@ public: return m_logA + m_b*logT - m_E*recipT; } + /** + * Update the value of the natural logarithm of the rate constant. + */ + doublereal updateLog(doublereal logT, doublereal recipT) const { + return m_logA + m_b*logT - m_E*recipT; + } + /** * Update the value the rate constant. * @@ -337,23 +344,13 @@ public: // 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]; - } + ihigh1_ = iter->second.first; + ihigh2_ = iter->second.second; // 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]; - } + ilow1_ = iter->second.first; + ilow2_ = iter->second.second; rDeltaP_ = 1.0 / (logP2_ - logP1_); } @@ -373,22 +370,22 @@ public: */ doublereal updateRC(doublereal logT, doublereal recipT) const { double log_k1, log_k2; - if (m1_ == 1) { - log_k1 = A1_[0] + n1_[0] * logT - Ea1_[0] * recipT; + if (ilow1_ == ilow2_) { + log_k1 = rates_[ilow1_].updateLog(logT, recipT); } else { double k = 1e-300; // non-zero to make log(k) finite - for (size_t m = 0; m < m1_; m++) { - k += A1_[m] * std::exp(n1_[m] * logT - Ea1_[m] * recipT); + for (size_t i = ilow1_; i < ilow2_; i++) { + k += rates_[i].updateRC(logT, recipT); } log_k1 = std::log(k); } - if (m2_ == 1) { - log_k2 = A2_[0] + n2_[0] * logT - Ea2_[0] * recipT; + if (ihigh1_ == ihigh2_) { + log_k2 = rates_[ihigh1_].updateLog(logT, recipT); } else { double k = 1e-300; // non-zero to make log(k) finite - for (size_t m = 0; m < m2_; m++) { - k += A2_[m] * std::exp(n2_[m] * logT - Ea2_[m] * recipT); + for (size_t i = ihigh1_; i < ihigh2_; i++) { + k += rates_[i].updateRC(logT, recipT); } log_k2 = std::log(k); } @@ -413,28 +410,23 @@ public: void validate(const std::string& equation); protected: - //! log(p) to (index range) in A_, n, Ea vectors + //! log(p) to (index range) in the rates_ vector 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] + // 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 - //! 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 + //! 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_; - //! 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 }; diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index f16967d21..c6fd84680 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -143,16 +143,13 @@ Plog::Plog(const ReactionData& rdata) : logP_(-1000) , logP1_(1000) , logP2_(-1000) - , m1_(npos) - , m2_(npos) , rDeltaP_(-1.0) - , maxRates_(1) { typedef std::multimap::const_iterator iter_t; size_t j = 0; - size_t rateCount = 0; // Insert intermediate pressures + rates_.reserve(rdata.plogParameters.size()); for (iter_t iter = rdata.plogParameters.begin(); iter != rdata.plogParameters.end(); iter++) { @@ -160,42 +157,20 @@ Plog::Plog(const ReactionData& rdata) 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] = std::log(A_[iter->second.first]); - } + rates_.push_back(Arrhenius(iter->second[0], iter->second[1], + iter->second[2])); } // 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_); - if (rdata.validate) { validate(rdata.equation); } @@ -205,13 +180,10 @@ Plog::Plog(const std::multimap& rates) : logP_(-1000) , logP1_(1000) , logP2_(-1000) - , m1_(npos) - , m2_(npos) , rDeltaP_(-1.0) - , maxRates_(1) { size_t j = 0; - size_t rateCount = 0; + rates_.reserve(rates.size()); // Insert intermediate pressures for (std::multimap::const_iterator iter = rates.begin(); iter != rates.end(); @@ -220,41 +192,18 @@ Plog::Plog(const std::multimap& rates) 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.preExponentialFactor()); - n_.push_back(iter->second.temperatureExponent()); - Ea_.push_back(iter->second.activationEnergy_R()); - } - - // 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] = std::log(A_[iter->second.first]); - } + rates_.push_back(iter->second); } // 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_); } void Plog::validate(const std::string& equation)