From 6f0fb5cf8df10905c852b825514ce2bee4dde6c8 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 29 Oct 2014 23:32:27 +0000 Subject: [PATCH] [Kinetics] Refactor and deprecate installReagents --- include/cantera/kinetics/AqueousKinetics.h | 2 - include/cantera/kinetics/GasKinetics.h | 2 - include/cantera/kinetics/InterfaceKinetics.h | 2 - include/cantera/kinetics/Kinetics.h | 2 + src/kinetics/AqueousKinetics.cpp | 67 +++----------- src/kinetics/GasKinetics.cpp | 66 +++----------- src/kinetics/InterfaceKinetics.cpp | 93 ++------------------ src/kinetics/Kinetics.cpp | 39 +++++++- 8 files changed, 67 insertions(+), 206 deletions(-) diff --git a/include/cantera/kinetics/AqueousKinetics.h b/include/cantera/kinetics/AqueousKinetics.h index 9a4835060..1245c5ec1 100644 --- a/include/cantera/kinetics/AqueousKinetics.h +++ b/include/cantera/kinetics/AqueousKinetics.h @@ -164,8 +164,6 @@ protected: private: void addElementaryReaction(ReactionData& r); - void installReagents(const ReactionData& r); - /** * Update the equilibrium constants in molar units. */ diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index 1910092bd..9ae3ec7ee 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -162,8 +162,6 @@ private: void addPlogReaction(ReactionData& r); void addChebyshevReaction(ReactionData& r); - virtual void installReagents(const ReactionData& r); - //! Update the equilibrium constants in molar units. void updateKc(); diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 5b0aab55a..0a511c2d1 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -307,8 +307,6 @@ public: void addGlobalReaction(ReactionData& r); - void installReagents(const ReactionData& r); - //! Update the equilibrium constants and stored electrochemical potentials //! in molar units for all reversible reactions and for all species. /*! diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index e000e5493..76ccfc136 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -779,6 +779,8 @@ public: */ virtual void addReaction(ReactionData& r); + //! @deprecated To be removed after Cantera 2.2. No longer called as part + //! of addReaction. virtual void installReagents(const ReactionData& r) { throw NotImplementedError("Kinetics::installReagents"); } diff --git a/src/kinetics/AqueousKinetics.cpp b/src/kinetics/AqueousKinetics.cpp index f90c3dcf6..4d86e6ce5 100644 --- a/src/kinetics/AqueousKinetics.cpp +++ b/src/kinetics/AqueousKinetics.cpp @@ -339,6 +339,16 @@ void AqueousKinetics::addReaction(ReactionData& r) if (r.reactionType == ELEMENTARY_RXN) { addElementaryReaction(r); } + + m_dn.push_back(accumulate(r.pstoich.begin(), r.pstoich.end(), 0.0) - + accumulate(r.rstoich.begin(), r.rstoich.end(), 0.0)); + if (r.reversible) { + m_revindex.push_back(nReactions()); + m_nrev++; + } else { + m_irrev.push_back(nReactions()); + m_nirrev++; + } Kinetics::addReaction(r); } @@ -348,63 +358,6 @@ void AqueousKinetics::addElementaryReaction(ReactionData& r) m_rates.install(nReactions(), r); } -void AqueousKinetics::installReagents(const ReactionData& r) -{ - size_t n, ns, m; - doublereal nsFlt; - doublereal reactantGlobalOrder = 0.0; - doublereal productGlobalOrder = 0.0; - size_t rnum = nReactions(); - - std::vector rk; - size_t nr = r.reactants.size(); - for (n = 0; n < nr; n++) { - nsFlt = r.rstoich[n]; - reactantGlobalOrder += nsFlt; - ns = (size_t) nsFlt; - if ((doublereal) ns != nsFlt) { - ns = std::max(ns, 1); - } - if (r.rstoich[n] != 0.0) { - m_rrxn[r.reactants[n]][rnum] += r.rstoich[n]; - } - for (m = 0; m < ns; m++) { - rk.push_back(r.reactants[n]); - } - } - m_reactants.push_back(rk); - - std::vector pk; - size_t np = r.products.size(); - for (n = 0; n < np; n++) { - nsFlt = r.pstoich[n]; - productGlobalOrder += nsFlt; - ns = (size_t) nsFlt; - if ((double) ns != nsFlt) { - ns = std::max(ns, 1); - } - if (r.pstoich[n] != 0.0) { - m_prxn[r.products[n]][rnum] += r.pstoich[n]; - } - for (m = 0; m < ns; m++) { - pk.push_back(r.products[n]); - } - } - m_products.push_back(pk); - - m_rxnstoich.add(nReactions(), r); - - if (r.reversible) { - m_dn.push_back(productGlobalOrder - reactantGlobalOrder); - m_revindex.push_back(nReactions()); - m_nrev++; - } else { - m_dn.push_back(productGlobalOrder - reactantGlobalOrder); - m_irrev.push_back(nReactions()); - m_nirrev++; - } -} - void AqueousKinetics::init() { m_kk = thermo().nSpecies(); diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index d700e77b6..a379161cd 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -451,6 +451,16 @@ void GasKinetics::addReaction(ReactionData& r) } // operations common to all reaction types + m_dn.push_back(accumulate(r.pstoich.begin(), r.pstoich.end(), 0.0) - + accumulate(r.rstoich.begin(), r.rstoich.end(), 0.0)); + + if (r.reversible) { + m_revindex.push_back(nReactions()); + m_nrev++; + } else { + m_irrev.push_back(nReactions()); + m_nirrev++; + } Kinetics::addReaction(r); } @@ -507,62 +517,6 @@ void GasKinetics::addChebyshevReaction(ReactionData& r) m_cheb_rates.install(nReactions(), r); } -void GasKinetics::installReagents(const ReactionData& r) -{ - size_t n, ns, m; - doublereal nsFlt; - doublereal reactantGlobalOrder = 0.0; - doublereal productGlobalOrder = 0.0; - size_t rnum = nReactions(); - - std::vector rk; - size_t nr = r.reactants.size(); - for (n = 0; n < nr; n++) { - nsFlt = r.rstoich[n]; - reactantGlobalOrder += nsFlt; - ns = (size_t) nsFlt; - if ((doublereal) ns != nsFlt) { - ns = std::max(ns, 1); - } - if (r.rstoich[n] != 0.0) { - m_rrxn[r.reactants[n]][rnum] += r.rstoich[n]; - } - for (m = 0; m < ns; m++) { - rk.push_back(r.reactants[n]); - } - } - m_reactants.push_back(rk); - - std::vector pk; - size_t np = r.products.size(); - for (n = 0; n < np; n++) { - nsFlt = r.pstoich[n]; - productGlobalOrder += nsFlt; - ns = (size_t) nsFlt; - if ((double) ns != nsFlt) { - ns = std::max(ns, 1); - } - if (r.pstoich[n] != 0.0) { - m_prxn[r.products[n]][rnum] += r.pstoich[n]; - } - for (m = 0; m < ns; m++) { - pk.push_back(r.products[n]); - } - } - m_products.push_back(pk); - m_rxnstoich.add(nReactions(), r); - - if (r.reversible) { - m_dn.push_back(productGlobalOrder - reactantGlobalOrder); - m_revindex.push_back(nReactions()); - m_nrev++; - } else { - m_dn.push_back(productGlobalOrder - reactantGlobalOrder); - m_irrev.push_back(nReactions()); - m_nirrev++; - } -} - void GasKinetics::init() { m_kk = thermo().nSpecies(); diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 60c6481dd..6790bdd87 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -798,6 +798,13 @@ void InterfaceKinetics::addReaction(ReactionData& r) addElementaryReaction(r); } + if (r.reversible) { + m_revindex.push_back(nReactions()); + m_nrev++; + } else { + m_irrev.push_back(nReactions()); + m_nirrev++; + } Kinetics::addReaction(r); m_rxnPhaseIsReactant.push_back(std::vector(nPhases(), false)); @@ -975,92 +982,6 @@ void InterfaceKinetics::setIOFlag(int ioFlag) } } -void InterfaceKinetics::installReagents(const ReactionData& r) -{ - size_t n, ns, m; - doublereal nsFlt; - - /* - * Obtain the current reaction index for the reaction that we - * are adding. The first reaction is labeled 0. - */ - size_t rnum = nReactions(); - - // vectors rk and pk are lists of species numbers, with - // repeated entries for species with stoichiometric - // coefficients > 1. This allows the reaction to be defined - // with unity reaction order for each reactant, and so the - // faster method 'multiply' can be used to compute the rate of - // progress instead of 'power'. - - std::vector rk; - size_t nr = r.reactants.size(); - for (n = 0; n < nr; n++) { - nsFlt = r.rstoich[n]; - ns = (size_t) nsFlt; - if ((doublereal) ns != nsFlt) { - ns = std::max(ns, 1); - } - /* - * Add to m_rrxn. m_rrxn is a vector of maps. m_rrxn has a length - * equal to the total number of species for each species, there - * exists a map, with the reaction number being the key, and the - * reactant stoichiometric coefficient being the value. - */ - m_rrxn[r.reactants[n]][rnum] = nsFlt; - for (m = 0; m < ns; m++) { - rk.push_back(r.reactants[n]); - } - } - /* - * Now that we have rk[], we add it into the vector m_reactants - * in the rnum index spot. Thus m_reactants[rnum] yields a vector - * of reactants for the rnum'th reaction - */ - m_reactants.push_back(rk); - std::vector pk; - size_t np = r.products.size(); - for (n = 0; n < np; n++) { - nsFlt = r.pstoich[n]; - ns = (size_t) nsFlt; - if ((doublereal) ns != nsFlt) { - ns = std::max(ns, 1); - } - /* - * Add to m_prxn. m_prxn is a vector of maps. m_prxn has a length - * equal to the total number of species for each species, there - * exists a map, with the reaction number being the key, and the - * product stoichiometric coefficient being the value. - */ - m_prxn[r.products[n]][rnum] = nsFlt; - for (m = 0; m < ns; m++) { - pk.push_back(r.products[n]); - } - } - /* - * Now that we have pk[], we add it into the vector m_products - * in the rnum index spot. Thus m_products[rnum] yields a vector - * of products for the rnum'th reaction - */ - m_products.push_back(pk); - /* - * Add this reaction to the stoichiometric coefficient manager. This - * calculates rates of species production from reaction rates of - * progress. - */ - m_rxnstoich.add(nReactions(), r); - /* - * register reaction in lists of reversible and irreversible rxns. - */ - if (r.reversible) { - m_revindex.push_back(nReactions()); - m_nrev++; - } else { - m_irrev.push_back(nReactions()); - m_nirrev++; - } -} - void InterfaceKinetics::addPhase(thermo_t& thermo) { Kinetics::addPhase(thermo); diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 1d1478b27..29c6e2339 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -334,7 +334,44 @@ void Kinetics::finalize() } void Kinetics::addReaction(ReactionData& r) { - installReagents(r); + // vectors rk and pk are lists of species numbers, with repeated entries + // for species with stoichiometric coefficients > 1. This allows the + // reaction to be defined with unity reaction order for each reactant, and + // so the faster method 'multiply' can be used to compute the rate of + // progress instead of 'power'. + std::vector rk; + 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) { + ns = std::max(ns, 1); + } + if (r.rstoich[n] != 0.0) { + m_rrxn[r.reactants[n]][m_ii] += r.rstoich[n]; + } + for (size_t m = 0; m < ns; m++) { + rk.push_back(r.reactants[n]); + } + } + m_reactants.push_back(rk); + + std::vector pk; + 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) { + ns = std::max(ns, 1); + } + if (r.pstoich[n] != 0.0) { + m_prxn[r.products[n]][m_ii] += r.pstoich[n]; + } + for (size_t m = 0; m < ns; m++) { + pk.push_back(r.products[n]); + } + } + m_products.push_back(pk); + m_rxnstoich.add(nReactions(), r); + installGroups(nReactions(), r.rgroups, r.pgroups); incrementRxnCount(); m_rxneqn.push_back(r.equation);