research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Kinetics] Simplify Plog by using Arrhenius objects internally

Browse source
This commit is contained in:
Ray Speth 2015-04-24 14:49:42 -04:00
parent e7946478df
commit 4151c10884
2 changed files with 32 additions and 91 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

62
include/cantera/kinetics/RxnRates.h
View file

@ -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<double, std::pair<size_t, size_t> > pressures_;
typedef std::map<double, std::pair<size_t, size_t> >::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<Arrhenius> 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
};

61
src/kinetics/RxnRates.cpp
View file

@ -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<double, vector_fp>::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<double, Arrhenius>& 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<double, Arrhenius>::const_iterator iter = rates.begin();
iter != rates.end();
@ -220,41 +192,18 @@ Plog::Plog(const std::multimap<double, Arrhenius>& 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)