[Kinetics] Prevent double counting in reaction path diagrams

This fixes the double counting that occurs in reactions like:

    H + HO2 => 2 OH

Fixes #377
This commit is contained in:
Ray Speth 2017-01-27 18:15:34 -05:00
parent 3093e6e6d4
commit 5a0fb579a8
2 changed files with 19 additions and 14 deletions

View file

@ -318,7 +318,12 @@ protected:
std::vector<vector_int> m_groups;
std::vector<Group> m_sgroup;
std::vector<std::string> m_elementSymbols;
//! m_transfer[reaction][reactant number][product number] where "reactant
//! number" means the number of the reactant in the reaction equation, e.g.
//! for "A+B -> C+D", "B" is reactant number 1 and "C" is product number 0.
std::map<size_t, std::map<size_t, std::map<size_t, Group> > > m_transfer;
std::vector<bool> m_determinate;
Array2D m_atoms;
std::map<std::string, size_t> m_enamemap;

View file

@ -444,17 +444,17 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s)
group_a0 = &r0;
group_b0 = &b0;
group_c0 = &p1;
m_transfer[i][kr0][kp0] = r0;
m_transfer[i][kr1][kp0] = b0;
m_transfer[i][kr1][kp1] = p1;
m_transfer[i][0][0] = r0;
m_transfer[i][1][0] = b0;
m_transfer[i][1][1] = p1;
} else {
group_a0 = &r1;
group_c0 = &p0;
b0 *= -1;
group_b0 = &b0;
m_transfer[i][kr1][kp1] = r1;
m_transfer[i][kr0][kp1] = b0;
m_transfer[i][kr0][kp0] = p0;
m_transfer[i][1][1] = r1;
m_transfer[i][0][1] = b0;
m_transfer[i][0][0] = p0;
}
logfile << " ";
group_a0->fmt(logfile, m_elementSymbols);
@ -479,9 +479,9 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s)
group_b1 = &b1;
group_c1 = &p0;
if (!b0.valid()) {
m_transfer[i][kr0][kp1] = r0;
m_transfer[i][kr1][kp1] = b0;
m_transfer[i][kr1][kp0] = p0;
m_transfer[i][0][1] = r0;
m_transfer[i][1][1] = b0;
m_transfer[i][1][0] = p0;
}
} else {
group_a1 = &r1;
@ -489,9 +489,9 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s)
b1 *= -1;
group_b1 = &b1;
if (!b0.valid()) {
m_transfer[i][kr1][kp0] = r1;
m_transfer[i][kr0][kp0] = b0;
m_transfer[i][kr0][kp1] = p1;
m_transfer[i][1][0] = r1;
m_transfer[i][0][0] = b0;
m_transfer[i][0][1] = p1;
}
}
logfile << " ";
@ -768,10 +768,10 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element,
}
f = 0.0;
} else {
if (!g[kkr][kkp]) {
if (!g[kr][kp]) {
f = 0.0;
} else {
f = g[kkr][kkp].nAtoms(m);
f = g[kr][kp].nAtoms(m);
}
}
} else {