From e2f15e97d457a9dbee690c541665f412f2e706b5 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 15 Nov 2014 00:46:54 +0000 Subject: [PATCH] [Kinetics] Add duplicate reaction checks using Reaction objects --- include/cantera/kinetics/Kinetics.h | 25 ++++++ src/kinetics/Kinetics.cpp | 127 ++++++++++++++++++++++++++++ src/kinetics/importKinetics.cpp | 3 + 3 files changed, 155 insertions(+) diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 5b2b39786..93e34f021 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -871,6 +871,15 @@ public: return false; } + //! Check for duplicate reactions. + /** + * If `throw_err` is true, then an exception will be thrown if any unmarked + * duplicate reactions are found. Otherwise, the indices of the first pair + * of duplicate reactions found will be returned. If no duplicate reactions + * are found, returns `(npos, npos)`. + */ + virtual std::pair checkDuplicates(bool throw_err=true) const; + /*! * Takes as input an array of properties for all species in the mechanism * and copies those values belonging to a particular phase to the output @@ -890,6 +899,22 @@ protected: throw NotImplementedError("Kinetics::updateROP"); } + //! Check whether `r1` and `r2` represent duplicate stoichiometries + //! This function returns a ratio if two reactions are duplicates of + //! one another, and 0.0 otherwise. + /*! + * `r1` and `r2` are maps of species key to stoichiometric coefficient, one + * for each reaction, where the species key is `1+k` for reactants and + * `-1-k` for products and `k` is the species index. + * + * @return 0.0 if the stoichiometries are not multiples of one another + * Otherwise, it returns the ratio of the stoichiometric coefficients. + * + * @ingroup kineticsmgr + */ + double checkDuplicateStoich(std::map& r1, + std::map& r2) const; + //! @name Stoichiometry management /*! * These objects and functions handle turning reaction extents into species diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index e8330ade9..df9baac8c 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -162,6 +162,133 @@ void Kinetics::assignShallowPointers(const std::vector & tpVector) } +std::pair Kinetics::checkDuplicates(bool throw_err) const +{ + //! Map of (key indicating participating species) to reaction numbers + std::map > participants; + std::vector > net_stoich; + + for (size_t i = 0; i < m_reactions.size(); i++) { + // Get data about this reaction + unsigned long int key = 0; + Reaction& R = *m_reactions[i]; + net_stoich.push_back(std::map()); + std::map& net = net_stoich.back(); + for (Composition::const_iterator iter = R.reactants.begin(); + iter != R.reactants.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + key += k*(k+1); + net[-1 -k] -= iter->second; + } + for (Composition::const_iterator iter = R.products.begin(); + iter != R.products.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + key += k*(k+1); + net[1+k] += iter->second; + } + + // Compare this reaction to others with similar participants + vector& related = participants[key]; + + for (size_t m = 0; m < related.size(); m++) { + Reaction& other = *m_reactions[related[m]]; + if (R.reaction_type != other.reaction_type) { + continue; // different reaction types + } else if (R.duplicate && other.duplicate) { + continue; // marked duplicates + } + doublereal c = checkDuplicateStoich(net_stoich[i], net_stoich[m]); + if (c == 0) { + continue; // stoichiometries differ (not by a multiple) + } else if (c < 0.0 && !R.reversible && !other.reversible) { + continue; // irreversible reactions in opposite directions + } else if (R.reaction_type == FALLOFF_RXN || + R.reaction_type == CHEMACT_RXN) { + ThirdBody& tb1 = dynamic_cast(R).third_body; + ThirdBody& tb2 = dynamic_cast(other).third_body; + bool thirdBodyOk = true; + for (size_t k = 0; k < nTotalSpecies(); k++) { + string s = kineticsSpeciesName(k); + if (tb1.efficiency(s) * tb2.efficiency(s) != 0.0) { + thirdBodyOk = false; + break; + } + } + if (thirdBodyOk) { + continue; // No overlap in third body efficiencies + } + } else if (R.reaction_type == THREE_BODY_RXN) { + ThirdBody& tb1 = dynamic_cast(R).third_body; + ThirdBody& tb2 = dynamic_cast(other).third_body; + bool thirdBodyOk = true; + for (size_t k = 0; k < nTotalSpecies(); k++) { + string s = kineticsSpeciesName(k); + if (tb1.efficiency(s) * tb2.efficiency(s) != 0.0) { + thirdBodyOk = false; + break; + } + } + if (thirdBodyOk) { + continue; // No overlap in third body efficiencies + } + } + if (throw_err) { + string msg = string("Undeclared duplicate reactions detected:\n") + +"Reaction "+int2str(i+1)+": "+other.equation() + +"\nReaction "+int2str(m+1)+": "+R.equation()+"\n"; + throw CanteraError("installReaction", msg); + } else { + return make_pair(i,m); + } + } + participants[key].push_back(i); + } + return make_pair(npos, npos); +} + +double Kinetics::checkDuplicateStoich(std::map& r1, + std::map& r2) const +{ + map::const_iterator b = r1.begin(), e = r1.end(); + int k1 = b->first; + // check for duplicate written in the same direction + doublereal ratio = 0.0; + if (r1[k1] && r2[k1]) { + ratio = r2[k1]/r1[k1]; + ++b; + bool different = false; + for (; b != e; ++b) { + k1 = b->first; + if (!r1[k1] || !r2[k1] || fabs(r2[k1]/r1[k1] - ratio) > 1.e-8) { + different = true; + break; + } + } + if (!different) { + return ratio; + } + } + + // check for duplicate written in the reverse direction + b = r1.begin(); + k1 = b->first; + if (r1[k1] == 0.0 || r2[-k1] == 0.0) { + return 0.0; + } + ratio = r2[-k1]/r1[k1]; + ++b; + for (; b != e; ++b) { + k1 = b->first; + if (!r1[k1] || !r2[-k1] || fabs(r2[-k1]/r1[k1] - ratio) > 1.e-8) { + return 0.0; + } + } + return ratio; +} + + void Kinetics::selectPhase(const doublereal* data, const thermo_t* phase, doublereal* phase_data) { diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index 46ba2c7cc..e406b275f 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -1188,6 +1188,9 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin, } } + if (check_for_duplicates) { + kin.checkDuplicates(); + } /* * Finalize the installation of the kinetics, now that we know * the true number of reactions in the mechanism, itot.