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

[Kinetics] Handle stoichiometry and reaction order in StoichManagerN

Browse source 
In the previous implementation, the higher-level classes were making unnecessary
assumptions about the optimizations used in StoichManagerN for reactions with
unity stoichometric coefficients.
This commit is contained in:
Ray Speth 2014-11-01 00:13:49 +00:00
parent 557ffffc5b
commit b2afb1a7a9
2 changed files with 23 additions and 44 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
include/cantera/kinetics/StoichManager.h
View file

@ -858,27 +858,37 @@ public:
}
bool frac = false;
for (size_t n = 0; n < stoich.size(); n++) {
if (stoich[n] != 1.0 || order[n] != 1.0) {
if (fmod(stoich[n], 1.0) || fmod(order[n], 1.0)) {
frac = true;
break;
}
}
if (frac) {
if (frac || k.size() > 3) {
m_loc[rxn] = m_cn_list.size();
m_cn_list.push_back(C_AnyN(rxn, k, order, stoich));
} else {
switch (k.size()) {
// Try to express the reaction with unity stoichiometric
// coefficients (by repeating species when necessary) so that the
// simpler 'multiply' function can be used to compute the rate
// instead of 'power'.
std::vector<size_t> kRep;
for (size_t n = 0; n < k.size(); n++) {
for (size_t i = 0; i < stoich[n]; i++)
kRep.push_back(k[n]);
}
switch (kRep.size()) {
case 1:
m_loc[rxn] = m_c1_list.size();
m_c1_list.push_back(C1(rxn, k[0]));
m_c1_list.push_back(C1(rxn, kRep[0]));
break;
case 2:
m_loc[rxn] = m_c2_list.size();
m_c2_list.push_back(C2(rxn, k[0], k[1]));
m_c2_list.push_back(C2(rxn, kRep[0], kRep[1]));
break;
case 3:
m_loc[rxn] = m_c3_list.size();
m_c3_list.push_back(C3(rxn, k[0], k[1], k[2]));
m_c3_list.push_back(C3(rxn, kRep[0], kRep[1], kRep[2]));
break;
default:
m_loc[rxn] = m_cn_list.size();

45
src/kinetics/Kinetics.cpp
View file

@ -379,12 +379,10 @@ void Kinetics::addReaction(ReactionData& r) {
// so the faster method 'multiply' can be used to compute the rate of
// progress instead of 'power'.
std::vector<size_t> rk;
bool fracReactants = false;
for (size_t n = 0; n < r.reactants.size(); n++) {
double nsFlt = r.rstoich[n];
size_t ns = (size_t) nsFlt;
if ((double) ns != nsFlt) {
fracReactants = true;
ns = std::max<size_t>(ns, 1);
}
if (r.rstoich[n] != 0.0) {
@ -397,12 +395,10 @@ void Kinetics::addReaction(ReactionData& r) {
m_reactants.push_back(rk);
std::vector<size_t> pk;
bool fracProducts = false;
for (size_t n = 0; n < r.products.size(); n++) {
double nsFlt = r.pstoich[n];
size_t ns = (size_t) nsFlt;
if ((double) ns != nsFlt) {
fracProducts = true;
ns = std::max<size_t>(ns, 1);
}
if (r.pstoich[n] != 0.0) {
@ -414,19 +410,14 @@ void Kinetics::addReaction(ReactionData& r) {
}
m_products.push_back(pk);
size_t irxn = nReactions();
bool doGlobal = false;
std::vector<size_t> extReactants = r.reactants;
vector_fp extRStoich = r.rstoich;
vector_fp extROrder = r.rorder;
// If we have a complete global reaction then we need to do something more
// complete than the previous treatment. Basically we will use the reactant
// manager to calculate the global forward reaction rate of progress.
// If the reaction order involves non-reactant species, add extra terms to
// the reactants with zero stoichiometry so that the stoichiometry manager
// can be used to compute the global forward reaction rate.
if (r.forwardFullOrder_.size() > 0) {
// Trigger a treatment where the order of the reaction and the
// stoichiometry are treated as different.
doGlobal = true;
size_t nsp = r.forwardFullOrder_.size();
// Set up a signal vector to indicate whether the species has been added
@ -459,34 +450,12 @@ void Kinetics::addReaction(ReactionData& r) {
}
}
// If the reaction is non-mass action add it in in a general way
// Reactants get extra terms for the forward reaction rate of progress
// that may have zero stoichiometries.
if (doGlobal) {
m_reactantStoich.add(irxn, extReactants, extROrder, extRStoich);
} else {
// this is confusing. The only issue should be whether rorder is different than rstoich!
if (fracReactants || r.global || rk.size() > 3) {
m_reactantStoich.add(irxn, r.reactants, r.rorder, r.rstoich);
} else {
m_reactantStoich.add(irxn, rk);
}
}
size_t irxn = nReactions();
m_reactantStoich.add(irxn, extReactants, extROrder, extRStoich);
if (r.reversible) {
// this is confusing. The only issue should be whether porder is different than pstoich!
if (pk.size() > 3 || r.isReversibleWithFrac) {
m_revProductStoich.add(irxn, r.products, r.porder, r.pstoich);
} else {
m_revProductStoich.add(irxn, pk);
}
m_revProductStoich.add(irxn, r.products, r.porder, r.pstoich);
} else {
// this is confusing. The only issue should be whether porder is different than pstoich!
if (fracProducts || pk.size() > 3) {
m_irrevProductStoich.add(irxn, r.products, r.porder, r.pstoich);
} else {
m_irrevProductStoich.add(irxn, pk);
}
m_irrevProductStoich.add(irxn, r.products, r.porder, r.pstoich);
}
installGroups(nReactions(), r.rgroups, r.pgroups);