[Kinetics] Enable calculations with ThreeBodyArrheniusRate

ischoegl · ischoegl · commit c93ea765fc37 · 2022-03-15T21:23:47.000-05:00
diff --git a/include/cantera/kinetics/BulkRate.h b/include/cantera/kinetics/BulkRate.h
@@ -69,6 +69,14 @@ class ThreeBodyBase
     //! Get the third-body efficiency for species *k*
     double efficiency(const std::string& k) const;
 
+    //! Get the third-body efficiency for species *k*
+    void getEfficiencyMap(std::map<size_t, double>& eff) const;
+
+    //! Get flag indicating whether third-body participates in the law of mass action
+    bool massAction() const {
+        return m_massAction;
+    }
+
     //! Build rate-specific parameters based on Reaction and Kinetics context
     //! @param rxn  Associated reaction object
     //! @param kin  Kinetics object holding the rate evaluator
@@ -79,7 +87,7 @@ class ThreeBodyBase
     void updateFromStruct(const BulkData& shared_data) {
         if (shared_data.ready) {
             m_thirdBodyConc = m_defaultEfficiency * shared_data.molarDensity;
-            for (const auto& eff : m_efficiencies) {
+            for (const auto& eff : m_efficiencyMap) {
                 m_thirdBodyConc += shared_data.concentrations[eff.first] * eff.second;
             }
         }
@@ -90,14 +98,6 @@ class ThreeBodyBase
         return m_thirdBodyConc;
     }
 
-    //! Apply correction
-    double applyCorrection(double value) const {
-        if (m_massAction) {
-            return value * m_thirdBodyConc;
-        }
-        return value;
-    }
-
 protected:
     double m_thirdBodyConc; //!< Effective third-body concentration
 
@@ -111,11 +111,11 @@ class ThreeBodyBase
     //! (`true` for three-body reactions, `false` for falloff reactions)
     bool m_massAction;
 
-    Composition m_efficiencyMap; //!< Map of species to third body efficiency
+    Composition m_efficiencies; //!< Composition defining third body efficiency
 
 private:
     //! Vector of pairs containing indices and efficiencies
-    std::vector<std::pair<size_t, double>> m_efficiencies;
+    std::vector<std::pair<size_t, double>> m_efficiencyMap;
 };
 
 
@@ -170,7 +170,7 @@ class ThreeBodyRate : public RateType, public ThreeBodyBase
     //! Evaluate reaction rate
     //! @param shared_data  data shared by all reactions of a given type
     double evalFromStruct(const DataType& shared_data) const {
-        return applyCorrection(RateType::evalFromStruct(shared_data));
+        return RateType::evalFromStruct(shared_data);
     }
 
 protected:
diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h
@@ -119,6 +119,10 @@ class Reaction
     //! @param detect  Run detection of unmarked third-body collider
     bool stripThirdBody(bool detect=true);
 
+    std::string thirdBodyCollider() const {
+        return m_thirdBodyCollider;
+    }
+
     //! Verify that all species involved in the reaction are defined in the Kinetics
     //! object. The function returns true if all species are found, and raises an
     //! exception unless the kinetics object is configured to skip undeclared species,
diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp
@@ -156,6 +156,15 @@ bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
         if (r->thirdBody() != nullptr) {
             addThirdBody(r);
         }
+        const auto threeBody = std::dynamic_pointer_cast<ThreeBodyArrheniusRate>(rate);
+        if (threeBody) {
+            // preliminary proof-of-concept uses ThirdBodyCalc3
+            std::map<size_t, double> efficiencies;
+            threeBody->getEfficiencyMap(efficiencies);
+            m_multi_concm.install(
+                nReactions() - 1, efficiencies,
+                threeBody->defaultEfficiency(), threeBody->massAction());
+        }
     }
 
     m_concm.push_back(NAN);
diff --git a/src/kinetics/BulkRate.cpp b/src/kinetics/BulkRate.cpp
@@ -21,16 +21,17 @@ void ThreeBodyBase::setParameters(const AnyMap& node)
 {
     m_defaultEfficiency = node.getDouble("default-efficiency", 1.0);
     if (node.hasKey("efficiencies")) {
-        m_efficiencyMap = node["efficiencies"].asMap<double>();
+        m_efficiencies = node["efficiencies"].asMap<double>();
     }
     m_specifiedCollisionPartner = node.getBool("specified-collider", false);
+    m_massAction = node.getBool("mass-action", true);
 }
 
 void ThreeBodyBase::getParameters(AnyMap& node) const
 {
     if (!m_specifiedCollisionPartner) {
-        if (m_efficiencyMap.size()) {
-            node["efficiencies"] = m_efficiencyMap;
+        if (m_efficiencies.size()) {
+            node["efficiencies"] = m_efficiencies;
             node["efficiencies"].setFlowStyle();
         }
         if (m_defaultEfficiency != 1.0) {
@@ -39,23 +40,32 @@ void ThreeBodyBase::getParameters(AnyMap& node) const
     }
 }
 
-void ThreeBodyBase::getEfficiencies(AnyMap& efficiencies) const {
+void ThreeBodyBase::getEfficiencies(AnyMap& efficiencies) const
+{
     efficiencies.clear();
-    for (const auto& eff : m_efficiencyMap) {
+    for (const auto& eff : m_efficiencies) {
         efficiencies[eff.first] = eff.second;
     }
 }
 
+void ThreeBodyBase::getEfficiencyMap(std::map<size_t, double>& eff) const
+{
+    eff.clear();
+    for (const auto& item : m_efficiencyMap) {
+        eff[item.first] = item.second + m_defaultEfficiency;
+    }
+}
+
 double ThreeBodyBase::efficiency(const std::string& k) const {
-    return getValue(m_efficiencyMap, k, m_defaultEfficiency);
+    return getValue(m_efficiencies, k, m_defaultEfficiency);
 }
 
 void ThreeBodyBase::setContext(const Reaction& rxn, const Kinetics& kin)
 {
-    for (const auto& eff : m_efficiencyMap) {
+    for (const auto& eff : m_efficiencies) {
         size_t k = kin.kineticsSpeciesIndex(eff.first);
         if (k != npos) {
-            m_efficiencies.emplace_back(k, eff.second - m_defaultEfficiency);
+            m_efficiencyMap.emplace_back(k, eff.second - m_defaultEfficiency);
         } else if (!kin.skipUndeclaredThirdBodies()) {
             throw CanteraError("ThreeBodyBase::setContext",
                 "Found third-body efficiency for undeclared species '{}' "

Original file line number	Diff line number	Diff line change
`@@ -21,16 +21,17 @@ void ThreeBodyBase::setParameters(const AnyMap& node)`
`21`	`21`	`{`
`22`	`22`	`m_defaultEfficiency = node.getDouble("default-efficiency", 1.0);`
`23`	`23`	`if (node.hasKey("efficiencies")) {`
`24`		`- m_efficiencyMap = node["efficiencies"].asMap<double>();`
	`24`	`+ m_efficiencies = node["efficiencies"].asMap<double>();`
`25`	`25`	`}`
`26`	`26`	`m_specifiedCollisionPartner = node.getBool("specified-collider", false);`
	`27`	`+ m_massAction = node.getBool("mass-action", true);`
`27`	`28`	`}`
`28`	`29`
`29`	`30`	`void ThreeBodyBase::getParameters(AnyMap& node) const`
`30`	`31`	`{`
`31`	`32`	`if (!m_specifiedCollisionPartner) {`
`32`		`- if (m_efficiencyMap.size()) {`
`33`		`- node["efficiencies"] = m_efficiencyMap;`
	`33`	`+ if (m_efficiencies.size()) {`
	`34`	`+ node["efficiencies"] = m_efficiencies;`
`34`	`35`	`node["efficiencies"].setFlowStyle();`
`35`	`36`	`}`
`36`	`37`	`if (m_defaultEfficiency != 1.0) {`
`@@ -39,23 +40,32 @@ void ThreeBodyBase::getParameters(AnyMap& node) const`
`39`	`40`	`}`
`40`	`41`	`}`
`41`	`42`
`42`		`-void ThreeBodyBase::getEfficiencies(AnyMap& efficiencies) const {`
	`43`	`+void ThreeBodyBase::getEfficiencies(AnyMap& efficiencies) const`
	`44`	`+{`
`43`	`45`	`efficiencies.clear();`
`44`		`- for (const auto& eff : m_efficiencyMap) {`
	`46`	`+ for (const auto& eff : m_efficiencies) {`
`45`	`47`	`efficiencies[eff.first] = eff.second;`
`46`	`48`	`}`
`47`	`49`	`}`
`48`	`50`
	`51`	`+void ThreeBodyBase::getEfficiencyMap(std::map<size_t, double>& eff) const`
	`52`	`+{`
	`53`	`+ eff.clear();`
	`54`	`+ for (const auto& item : m_efficiencyMap) {`
	`55`	`+ eff[item.first] = item.second + m_defaultEfficiency;`
	`56`	`+ }`
	`57`	`+}`
	`58`	`+`
`49`	`59`	`double ThreeBodyBase::efficiency(const std::string& k) const {`
`50`		`- return getValue(m_efficiencyMap, k, m_defaultEfficiency);`
	`60`	`+ return getValue(m_efficiencies, k, m_defaultEfficiency);`
`51`	`61`	`}`
`52`	`62`
`53`	`63`	`void ThreeBodyBase::setContext(const Reaction& rxn, const Kinetics& kin)`
`54`	`64`	`{`
`55`		`- for (const auto& eff : m_efficiencyMap) {`
	`65`	`+ for (const auto& eff : m_efficiencies) {`
`56`	`66`	`size_t k = kin.kineticsSpeciesIndex(eff.first);`
`57`	`67`	`if (k != npos) {`
`58`		`- m_efficiencies.emplace_back(k, eff.second - m_defaultEfficiency);`
	`68`	`+ m_efficiencyMap.emplace_back(k, eff.second - m_defaultEfficiency);`
`59`	`69`	`} else if (!kin.skipUndeclaredThirdBodies()) {`
`60`	`70`	`throw CanteraError("ThreeBodyBase::setContext",`
`61`	`71`	`"Found third-body efficiency for undeclared species '{}' "`