From 62cd3250365efaa59fa1c8e3f878827ae360be75 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 7 May 2015 18:31:40 -0400 Subject: [PATCH] [Kinetics] Implement modifyReaction for gas-phase reactions --- include/cantera/kinetics/BulkKinetics.h | 1 + include/cantera/kinetics/FalloffMgr.h | 14 +++ include/cantera/kinetics/GasKinetics.h | 11 ++ include/cantera/kinetics/Kinetics.h | 10 ++ include/cantera/kinetics/RateCoeffMgr.h | 10 ++ interfaces/cython/cantera/_cantera.pxd | 1 + interfaces/cython/cantera/kinetics.pyx | 15 +++ .../cython/cantera/test/test_kinetics.py | 103 ++++++++++++++++++ src/kinetics/BulkKinetics.cpp | 5 + src/kinetics/GasKinetics.cpp | 58 ++++++++++ src/kinetics/Kinetics.cpp | 24 ++++ 11 files changed, 252 insertions(+) diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h index 720327a4b..c3d1c9795 100644 --- a/include/cantera/kinetics/BulkKinetics.h +++ b/include/cantera/kinetics/BulkKinetics.h @@ -46,6 +46,7 @@ public: protected: virtual void addElementaryReaction(ReactionData& r); virtual void addElementaryReaction(ElementaryReaction& r); + virtual void modifyElementaryReaction(size_t i, ElementaryReaction& rNew); Rate1 m_rates; std::vector m_revindex; //!< Indices of reversible reactions diff --git a/include/cantera/kinetics/FalloffMgr.h b/include/cantera/kinetics/FalloffMgr.h index f3a98423a..a1804a9da 100644 --- a/include/cantera/kinetics/FalloffMgr.h +++ b/include/cantera/kinetics/FalloffMgr.h @@ -61,6 +61,17 @@ public: m_worksize += f->workSize(); m_falloff.push_back(f); m_reactionType.push_back(reactionType); + m_indices[rxn] = m_falloff.size()-1; + } + + /*! + * Replace an existing falloff function calculator + * + * @param rxn External reaction index + * @param f New falloff function, of the same kind as the existing one + */ + void replace(size_t rxn, shared_ptr f) { + m_falloff[m_indices[rxn]] = f; } //! Size of the work array required to store intermediate results. @@ -109,6 +120,9 @@ protected: //! Distinguish between falloff and chemically activated reactions vector_int m_reactionType; + + //! map of external reaction index to local index + std::map m_indices; }; } diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index de39b2806..5951b0bed 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -54,6 +54,7 @@ public: virtual void init(); virtual void addReaction(ReactionData& r); virtual bool addReaction(shared_ptr r); + virtual void modifyReaction(size_t i, shared_ptr rNew); virtual void finalize(); virtual bool ready() const; //@} @@ -73,8 +74,13 @@ public: protected: size_t m_nfall; + //! Reaction index of each falloff reaction std::vector m_fallindx; + //! Map of reaction index to falloff reaction index (i.e indices in + //! #m_falloff_low_rates and #m_falloff_high_rates) + std::map m_rfallindx; + Rate1 m_falloff_low_rates; Rate1 m_falloff_high_rates; @@ -112,6 +118,11 @@ protected: void addPlogReaction(PlogReaction& r); void addChebyshevReaction(ChebyshevReaction& r); + void modifyThreeBodyReaction(size_t i, ThirdBodyReaction& r); + void modifyFalloffReaction(size_t i, FalloffReaction& r); + void modifyPlogReaction(size_t i, PlogReaction& r); + void modifyChebyshevReaction(size_t i, ChebyshevReaction& r); + //! Update the equilibrium constants in molar units. void updateKc(); diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 40748a609..5992579e1 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -807,6 +807,16 @@ public: */ virtual bool addReaction(shared_ptr r); + /** + * Modify the rate expression associated with a reaction. The + * stoichiometric equation, type of the reaction, reaction orders, third + * body efficiencies, reversibility, etc. must be unchanged. + * + * @param i Index of the reaction to be modified + * @param rNew Reaction with the new rate expressions + */ + virtual void modifyReaction(size_t i, shared_ptr rNew); + /** * Return the Reaction object for reaction *i*. */ diff --git a/include/cantera/kinetics/RateCoeffMgr.h b/include/cantera/kinetics/RateCoeffMgr.h index 2898ea98c..f477cba06 100644 --- a/include/cantera/kinetics/RateCoeffMgr.h +++ b/include/cantera/kinetics/RateCoeffMgr.h @@ -42,6 +42,7 @@ public: // Install a rate calculator and return the index of the calculator. m_rxn.push_back(rxnNumber); m_rates.push_back(R(rdata)); + m_indices[rxnNumber] = m_rxn.size() - 1; return m_rates.size() - 1; } @@ -53,8 +54,14 @@ public: void install(size_t rxnNumber, const R& rate) { m_rxn.push_back(rxnNumber); m_rates.push_back(rate); + m_indices[rxnNumber] = m_rxn.size() - 1; } + //! Replace an existing rate coefficient calculator + void replace(size_t rxnNumber, const R& rate) { + size_t i = m_indices[rxnNumber]; + m_rates[i] = rate; + } /** * Update the concentration-dependent parts of the rate @@ -94,6 +101,9 @@ public: protected: std::vector m_rates; std::vector m_rxn; + + //! map reaction number to index in m_rxn / m_rates + std::map m_indices; }; } diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 92e2a5dc3..9daf99dd7 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -342,6 +342,7 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera": void skipUndeclaredThirdBodies(cbool) void addReaction(shared_ptr[CxxReaction]) except + void finalize() except + + void modifyReaction(int, shared_ptr[CxxReaction]) except + shared_ptr[CxxReaction] reaction(size_t) except + cbool isReversible(int) except + diff --git a/interfaces/cython/cantera/kinetics.pyx b/interfaces/cython/cantera/kinetics.pyx index 70f648d5a..e7956c73e 100644 --- a/interfaces/cython/cantera/kinetics.pyx +++ b/interfaces/cython/cantera/kinetics.pyx @@ -74,8 +74,23 @@ cdef class Kinetics(_SolutionBase): return self.kinetics.kineticsSpeciesIndex(k, phase) def reaction(self, int i_reaction): + """ + Return a `Reaction` object representing the reaction with index + ``i_reaction``. + """ return wrapReaction(self.kinetics.reaction(i_reaction)) + def modify_reaction(self, int irxn, Reaction rxn): + """ + Modify the `Reaction` with index ``irxn`` to have the same rate + parameters as ``rxn``. ``rxn`` must have the same reactants and products + and be of the same type (i.e. `ElementaryReaction`, `FalloffReaction`, + `PlogReaction`, etc.) as the existing reaction. This method does not + modify the third-body efficiencies, reaction orders, or reversibility of + the reaction. + """ + self.kinetics.modifyReaction(irxn, rxn._reaction) + def is_reversible(self, int i_reaction): """True if reaction `i_reaction` is reversible.""" self._check_reaction_index(i_reaction) diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 164527a02..b57c5a9e1 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -765,3 +765,106 @@ class TestReaction(utilities.CanteraTest): surf2.forward_rate_constants[0]) self.assertNear(surf1.net_rates_of_progress[1], surf2.net_rates_of_progress[0]) + + def test_modify_invalid(self): + # different reaction type + tbr = self.gas.reaction(0) + R2 = ct.ElementaryReaction(tbr.reactants, tbr.products) + R2.rate = tbr.rate + with self.assertRaises(Exception): + self.gas.modify_reaction(0, R2) + + # different reactants + R = self.gas.reaction(7) + with self.assertRaises(Exception): + self.gas.modify_reaction(23, R) + + # different products + R = self.gas.reaction(14) + with self.assertRaises(Exception): + self.gas.modify_reaction(15, R) + + def test_modify_elementary(self): + gas = ct.Solution('h2o2.xml') + gas.TPX = self.gas.TPX + R = self.gas.reaction(2) + A1 = R.rate.preexponential_factor + b1 = R.rate.temperature_exponent + Ta1 = R.rate.activation_energy / ct.gas_constant + T = gas.T + self.assertNear(A1*T**b1*np.exp(-Ta1/T), gas.forward_rate_constants[2]) + + A2 = 1.5 * A1 + b2 = b1 + 0.1 + Ta2 = Ta1 * 1.2 + R.rate = ct.Arrhenius(A2, b2, Ta2 * ct.gas_constant) + gas.modify_reaction(2, R) + self.assertNear(A2*T**b2*np.exp(-Ta2/T), gas.forward_rate_constants[2]) + + def test_modify_third_body(self): + gas = ct.Solution('h2o2.xml') + gas.TPX = self.gas.TPX + R = self.gas.reaction(5) + A1 = R.rate.preexponential_factor + b1 = R.rate.temperature_exponent + T = gas.T + kf1 = gas.forward_rate_constants[5] + + A2 = 1.7 * A1 + b2 = b1 - 0.1 + R.rate = ct.Arrhenius(A2, b2, 0.0) + gas.modify_reaction(5, R) + kf2 = gas.forward_rate_constants[5] + self.assertNear((A2*T**b2) / (A1*T**b1), kf2/kf1) + + def test_modify_falloff(self): + gas = ct.Solution('gri30.xml') + gas.TPX = 1100, 3 * ct.one_atm, 'CH4:1.0, O2:0.4, CO2:0.1, H2O:0.05' + # these two reactions happen to have the same third-body efficiencies + r1 = gas.reaction(49) + r2 = gas.reaction(53) + self.assertEqual(r1.efficiencies, r2.efficiencies) + r2.high_rate = r1.high_rate + r2.low_rate = r1.low_rate + r2.falloff = r1.falloff + + gas.modify_reaction(53, r2) + kf = gas.forward_rate_constants + self.assertNear(kf[49], kf[53]) + + def test_modify_plog(self): + gas = ct.Solution('pdep-test.cti') + gas.TPX = 1010, 0.12 * ct.one_atm, 'R1A:0.3, R1B:0.2, H:0.1, R2:0.4' + + r0 = gas.reaction(0) + r1 = gas.reaction(1) + r0.rates = r1.rates + gas.modify_reaction(0, r0) + kf = gas.forward_rate_constants + self.assertNear(kf[0], kf[1]) + + # Removing the high-pressure rates should have no effect at low P... + r1.rates = r1.rates[:-4] + gas.modify_reaction(1, r1) + self.assertNear(kf[1], gas.forward_rate_constants[1]) + + # ... but should change the rate at higher pressures + gas.TP = 1010, 12.0 * ct.one_atm + kf = gas.forward_rates_of_progress + self.assertNotAlmostEqual(kf[0], kf[1]) + + def test_modify_chebyshev(self): + gas = ct.Solution('pdep-test.cti') + gas.TPX = 1010, 0.34 * ct.one_atm, 'R1A:0.3, R1B:0.2, H:0.1, R2:0.4' + + r1 = gas.reaction(4) + r2 = gas.reaction(5) + r1.set_parameters(r2.Tmin, r2.Tmax, r2.Pmin, r2.Pmax, r2.coeffs) + + # rates should be different before calling 'modify_reaction' + kf = gas.forward_rate_constants + self.assertNotAlmostEqual(kf[4], kf[5]) + + gas.modify_reaction(4, r1) + kf = gas.forward_rate_constants + self.assertNear(kf[4], kf[5]) diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index 6f43be179..f7ba6a8f2 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -166,6 +166,11 @@ void BulkKinetics::addElementaryReaction(ElementaryReaction& r) m_rates.install(nReactions()-1, r.rate); } +void BulkKinetics::modifyElementaryReaction(size_t i, ElementaryReaction& rNew) +{ + m_rates.replace(i, rNew.rate); +} + void BulkKinetics::init() { m_kk = thermo().nSpecies(); diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 65aaad282..df4f8d6c5 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -311,6 +311,7 @@ void GasKinetics::addFalloffReaction(ReactionData& r) // add this reaction number to the list of falloff reactions m_fallindx.push_back(nReactions()); + m_rfallindx[nReactions()] = m_nfall; // install the enhanced third-body concentration calculator for this // reaction @@ -353,6 +354,7 @@ void GasKinetics::addFalloffReaction(FalloffReaction& r) // add this reaction number to the list of falloff reactions m_fallindx.push_back(nReactions()-1); + m_rfallindx[nReactions()-1] = m_nfall; // install the enhanced third-body concentration calculator map efficiencies; @@ -408,6 +410,62 @@ void GasKinetics::addChebyshevReaction(ChebyshevReaction& r) m_cheb_rates.install(nReactions()-1, r.rate); } +void GasKinetics::modifyReaction(size_t i, shared_ptr rNew) +{ + // operations common to all reaction types + BulkKinetics::modifyReaction(i, rNew); + + switch (rNew->reaction_type) { + case ELEMENTARY_RXN: + modifyElementaryReaction(i, dynamic_cast(*rNew)); + break; + case THREE_BODY_RXN: + modifyThreeBodyReaction(i, dynamic_cast(*rNew)); + break; + case FALLOFF_RXN: + case CHEMACT_RXN: + modifyFalloffReaction(i, dynamic_cast(*rNew)); + break; + case PLOG_RXN: + modifyPlogReaction(i, dynamic_cast(*rNew)); + break; + case CHEBYSHEV_RXN: + modifyChebyshevReaction(i, dynamic_cast(*rNew)); + break; + default: + throw CanteraError("GasKinetics::modifyReaction", + "Unknown reaction type specified: " + int2str(rNew->reaction_type)); + } + + // invalidate all cached data + m_ROP_ok = false; + m_temp += 0.1234; + m_pres += 0.1234; +} + +void GasKinetics::modifyThreeBodyReaction(size_t i, ThirdBodyReaction& r) +{ + m_rates.replace(i, r.rate); +} + +void GasKinetics::modifyFalloffReaction(size_t i, FalloffReaction& r) +{ + size_t iFall = m_rfallindx[i]; + m_falloff_high_rates.replace(iFall, r.high_rate); + m_falloff_low_rates.replace(iFall, r.low_rate); + m_falloffn.replace(iFall, r.falloff); +} + +void GasKinetics::modifyPlogReaction(size_t i, PlogReaction& r) +{ + m_plog_rates.replace(i, r.rate); +} + +void GasKinetics::modifyChebyshevReaction(size_t i, ChebyshevReaction& r) +{ + m_cheb_rates.replace(i, r.rate); +} + void GasKinetics::init() { BulkKinetics::init(); diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 7b4aeeb1c..4bbdb57ad 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -770,6 +770,30 @@ bool Kinetics::addReaction(shared_ptr r) return true; } +void Kinetics::modifyReaction(size_t i, shared_ptr rNew) +{ + checkReactionIndex(i); + shared_ptr& rOld = m_reactions[i]; + if (rNew->reaction_type != rOld->reaction_type) { + throw CanteraError("Kinetics::modifyReaction", + "Reaction types are different: " + int2str(rOld->reaction_type) + + " != " + int2str(rNew->reaction_type) + "."); + } + + if (rNew->reactants != rOld->reactants) { + throw CanteraError("Kinetics::modifyReaction", + "Reactants are different: '" + rOld->reactantString() + "' != '" + + rNew->reactantString() + "'."); + } + + if (rNew->products != rOld->products) { + throw CanteraError("Kinetics::modifyReaction", + "Products are different: '" + rOld->productString() + "' != '" + + rNew->productString() + "'."); + } + m_reactions[i] = rNew; +} + shared_ptr Kinetics::reaction(size_t i) { checkReactionIndex(i);