From b5582e849f5a9e0ce1ac800bebeb45eb25677ec2 Mon Sep 17 00:00:00 2001
From: Raymond Speth <speth@mit.edu>
Date: Wed, 27 Jan 2021 18:43:41 -0500
Subject: [PATCH] [Kinetics] Fix third body concentrations for non-ideal gases

From the discussion on #967, the calcuation of [M] should be based on
the physical concentrations of the 3rd body colliders, not their
"activity concentrations".
---
 include/cantera/kinetics/BulkKinetics.h |  7 ++++++-
 include/cantera/kinetics/GasKinetics.h  |  1 +
 src/kinetics/BulkKinetics.cpp           |  3 ++-
 src/kinetics/GasKinetics.cpp            | 11 ++++++-----
 4 files changed, 15 insertions(+), 7 deletions(-)

diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h
index a778e199e15..471f3bb5109 100644
--- a/include/cantera/kinetics/BulkKinetics.h
+++ b/include/cantera/kinetics/BulkKinetics.h
@@ -62,7 +62,12 @@ class BulkKinetics : public Kinetics
     //! valued stoichiometries.
     vector_fp m_dn;
 
-    vector_fp m_conc;
+    //! Activity concentrations, as calculated by
+    //! ThermoPhase::getActivityConcentrations
+    vector_fp m_act_conc;
+
+    //! Physical concentrations, as calculated by ThermoPhase::getConcentrations
+    vector_fp m_phys_conc;
     vector_fp m_grt;
 
     bool m_ROP_ok;
diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h
index 1973323a7c4..bb1bf0ecdb9 100644
--- a/include/cantera/kinetics/GasKinetics.h
+++ b/include/cantera/kinetics/GasKinetics.h
@@ -101,6 +101,7 @@ class GasKinetics : public BulkKinetics
     vector_fp falloff_work;
     vector_fp concm_3b_values;
     vector_fp concm_falloff_values;
+
     //!@}
 
     void processFalloffReactions();
diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp
index f85dcba8180..270545709b9 100644
--- a/src/kinetics/BulkKinetics.cpp
+++ b/src/kinetics/BulkKinetics.cpp
@@ -178,7 +178,8 @@ void BulkKinetics::modifyElementaryReaction(size_t i, ElementaryReaction& rNew)
 void BulkKinetics::resizeSpecies()
 {
     Kinetics::resizeSpecies();
-    m_conc.resize(m_kk);
+    m_act_conc.resize(m_kk);
+    m_phys_conc.resize(m_kk);
     m_grt.resize(m_kk);
 }
 
diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp
index ef008ab112f..b40f95431df 100644
--- a/src/kinetics/GasKinetics.cpp
+++ b/src/kinetics/GasKinetics.cpp
@@ -68,17 +68,18 @@ void GasKinetics::update_rates_T()
 
 void GasKinetics::update_rates_C()
 {
-    thermo().getActivityConcentrations(m_conc.data());
+    thermo().getActivityConcentrations(m_act_conc.data());
+    thermo().getConcentrations(m_phys_conc.data());
     doublereal ctot = thermo().molarDensity();
 
     // 3-body reactions
     if (!concm_3b_values.empty()) {
-        m_3b_concm.update(m_conc, ctot, concm_3b_values.data());
+        m_3b_concm.update(m_phys_conc, ctot, concm_3b_values.data());
     }
 
     // Falloff reactions
     if (!concm_falloff_values.empty()) {
-        m_falloff_concm.update(m_conc, ctot, concm_falloff_values.data());
+        m_falloff_concm.update(m_phys_conc, ctot, concm_falloff_values.data());
     }
 
     // P-log reactions
@@ -187,10 +188,10 @@ void GasKinetics::updateROP()
     }
 
     // multiply ropf by concentration products
-    m_reactantStoich.multiply(m_conc.data(), m_ropf.data());
+    m_reactantStoich.multiply(m_act_conc.data(), m_ropf.data());
 
     // for reversible reactions, multiply ropr by concentration products
-    m_revProductStoich.multiply(m_conc.data(), m_ropr.data());
+    m_revProductStoich.multiply(m_act_conc.data(), m_ropr.data());
 
     for (size_t j = 0; j != nReactions(); ++j) {
         m_ropnet[j] = m_ropf[j] - m_ropr[j];