[Kinetics] Implement modifyReaction for gas-phase reactions

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<Arrhenius> m_rates;
     std::vector<size_t> m_revindex; //!< Indices of reversible reactions

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<Falloff> 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<size_t, size_t> m_indices;
 };
 }
 

include/cantera/kinetics/GasKinetics.h

@@ -54,6 +54,7 @@ public:
     virtual void init();
     virtual void addReaction(ReactionData& r);
     virtual bool addReaction(shared_ptr<Reaction> r);
+    virtual void modifyReaction(size_t i, shared_ptr<Reaction> 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<size_t> 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<size_t, size_t> m_rfallindx;
+
     Rate1<Arrhenius>                    m_falloff_low_rates;
     Rate1<Arrhenius>                    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();
 

include/cantera/kinetics/Kinetics.h

@@ -807,6 +807,16 @@ public:
      */
     virtual bool addReaction(shared_ptr<Reaction> 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<Reaction> rNew);
+
     /**
      * Return the Reaction object for reaction *i*.
      */

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<R>             m_rates;
     std::vector<size_t>           m_rxn;
+
+    //! map reaction number to index in m_rxn / m_rates
+    std::map<size_t, size_t> m_indices;
 };
 
 }

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 +

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)

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])

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();

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<size_t, double> 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<Reaction> rNew)
+{
+    // operations common to all reaction types
+    BulkKinetics::modifyReaction(i, rNew);
+
+    switch (rNew->reaction_type) {
+    case ELEMENTARY_RXN:
+        modifyElementaryReaction(i, dynamic_cast<ElementaryReaction&>(*rNew));
+        break;
+    case THREE_BODY_RXN:
+        modifyThreeBodyReaction(i, dynamic_cast<ThirdBodyReaction&>(*rNew));
+        break;
+    case FALLOFF_RXN:
+    case CHEMACT_RXN:
+        modifyFalloffReaction(i, dynamic_cast<FalloffReaction&>(*rNew));
+        break;
+    case PLOG_RXN:
+        modifyPlogReaction(i, dynamic_cast<PlogReaction&>(*rNew));
+        break;
+    case CHEBYSHEV_RXN:
+        modifyChebyshevReaction(i, dynamic_cast<ChebyshevReaction&>(*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();

src/kinetics/Kinetics.cpp

@@ -770,6 +770,30 @@ bool Kinetics::addReaction(shared_ptr<Reaction> r)
     return true;
 }
 
+void Kinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
+{
+    checkReactionIndex(i);
+    shared_ptr<Reaction>& 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<Reaction> Kinetics::reaction(size_t i)
 {
     checkReactionIndex(i);