Modify third body efficiencies using modify_reaction

include/cantera/kinetics/ThirdBodyCalc.h

@@ -32,6 +32,22 @@ class ThirdBodyCalc
             m_eff.back().push_back(eff.second - dflt);
         }
     }
+
+    void modifyEfficiencies(size_t rxnNumber, const std::map<size_t, double>& enhanced,
+                 double dflt=1.0) {
+        auto it = find(m_reaction_index.begin(), m_reaction_index.end(), rxnNumber);
+        if (it != m_reaction_index.end()) {
+            int index = it - m_reaction_index.begin();  
+            for (const auto& eff : enhanced) {
+                auto si = find(m_species[index].begin(), 
+                               m_species[index].end(), eff.first);   
+                if (it != m_reaction_index.end()) {
+                    int j = si - m_species[index].begin();
+                    m_eff[index][j] = (eff.second - m_default[index]);
+                }
+            }
+        }
+    }
 
     void update(const vector_fp& conc, double ctot, double* work) {
         for (size_t i = 0; i < m_species.size(); i++) {

src/kinetics/GasKinetics.cpp

@@ -395,6 +395,16 @@ void GasKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
 void GasKinetics::modifyThreeBodyReaction(size_t i, ThreeBodyReaction2& r)
 {
     m_rates.replace(i, r.rate);
+
+    map<size_t, double> efficiencies;
+    for (const auto& eff : r.third_body.efficiencies) {
+        size_t k = kineticsSpeciesIndex(eff.first);
+        if (k != npos) {
+            efficiencies[k] = eff.second;
+        }
+    }
+    m_3b_concm.modifyEfficiencies(i, efficiencies,
+                       r.third_body.default_efficiency);
 }
 
 void GasKinetics::modifyFalloffReaction(size_t i, FalloffReaction& r)
@@ -403,6 +413,16 @@ void GasKinetics::modifyFalloffReaction(size_t i, FalloffReaction& r)
     m_falloff_high_rates.replace(iFall, r.high_rate);
     m_falloff_low_rates.replace(iFall, r.low_rate);
     m_falloffn.replace(iFall, r.falloff);
+
+    map<size_t, double> efficiencies;
+    for (const auto& eff : r.third_body.efficiencies) {
+        size_t k = kineticsSpeciesIndex(eff.first);
+        if (k != npos) {
+            efficiencies[k] = eff.second;
+        }
+    }
+    m_falloff_concm.modifyEfficiencies(iFall, efficiencies,
+                       r.third_body.default_efficiency);
 }
 
 void GasKinetics::modifyPlogReaction(size_t i, PlogReaction2& r)

test/kinetics/kineticsFromScratch.cpp

@@ -46,6 +46,24 @@ class KineticsFromScratch : public testing::Test
         kin_ref->getRevRateConstants(&k_ref[0]);
         EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]);
     }
+
+    void check_uneq_rates(int iRef) {
+        ASSERT_EQ((size_t) 1, kin.nReactions());
+
+        std::string X = "O:0.02 H2:0.2 O2:0.5 H:0.03 OH:0.05 H2O:0.1 HO2:0.01";
+        p.setState_TPX(1200, 5*OneAtm, X);
+        p_ref.setState_TPX(1200, 5*OneAtm, X);
+
+        vector_fp k(1), k_ref(kin_ref->nReactions());
+
+        kin.getFwdRateConstants(&k[0]);
+        kin_ref->getFwdRateConstants(&k_ref[0]);
+        ASSERT_NE(k_ref[iRef], k[0]);
+
+        kin.getRevRateConstants(&k[0]);
+        kin_ref->getRevRateConstants(&k_ref[0]);
+        ASSERT_NE(k_ref[iRef], k[0]);
+    }
 };
 
 TEST_F(KineticsFromScratch, add_elementary_reaction)
@@ -118,10 +136,18 @@ TEST_F(KineticsFromScratch, add_falloff_reaction)
     Arrhenius low_rate(2.3e12, -0.9, -1700.0 / GasConst_cal_mol_K);
     vector_fp falloff_params { 0.7346, 94.0, 1756.0, 5182.0 };
     ThirdBody tbody;
-    tbody.efficiencies = parseCompString("AR:0.7 H2:2.0 H2O:6.0");
+    tbody.efficiencies = parseCompString("AR:0.7 H2:2.0 H2O:60.0");
     auto R = make_shared<FalloffReaction>(reac, prod, low_rate, high_rate, tbody);
     R->falloff = newFalloff("Troe", falloff_params);
+
+    ThirdBody tbody_correct;
+    tbody_correct.efficiencies = parseCompString("AR:0.7 H2:2.0 H2O:6.0");
+    auto R_correct = make_shared<FalloffReaction>(reac, prod, low_rate, high_rate, tbody_correct);
+    R_correct->falloff = newFalloff("Troe", falloff_params);
+
     kin.addReaction(R);
+    check_uneq_rates(2);
+    kin.modifyReaction(0, R_correct);
     check_rates(2);
 }