From 23642ded3c0bae446b0644adf136b5f6e6ca4f0a Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Wed, 29 Apr 2015 00:29:40 -0400
Subject: [PATCH] [Kinetics] Use current site density sticking reaction rate
 coefficient

The conversion of sticking coefficients to Arrhenius expressions includes the
site density on the surface phase, which may not be set when the reaction is
being added to the Kinetics object. We now defer the calculation of this
dependency until the rate coefficient is being evaluated so that the current
value for the site density can be used.
---
 include/cantera/kinetics/InterfaceKinetics.h |  7 ++++
 include/cantera/kinetics/Kinetics.h          |  3 ++
 src/kinetics/InterfaceKinetics.cpp           | 40 +++++++++++++++++++-
 3 files changed, 48 insertions(+), 2 deletions(-)

diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h
index f4da0f3ac..21d2dce02 100644
--- a/include/cantera/kinetics/InterfaceKinetics.h
+++ b/include/cantera/kinetics/InterfaceKinetics.h
@@ -719,6 +719,13 @@ protected:
      */
     std::vector<std::vector<bool> > m_rxnPhaseIsProduct;
 
+    //! Pairs of (reaction index, total order) for sticking reactions, which are
+    //! needed to compute the dependency of the rate constant on the site
+    //! density.
+    std::vector<std::pair<size_t, double> > m_sticking_orders;
+
+    void applyStickingCorrection(double* kf);
+
     int m_ioFlag;
 };
 }
diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h
index 4ef547ad8..7d4174cb7 100644
--- a/include/cantera/kinetics/Kinetics.h
+++ b/include/cantera/kinetics/Kinetics.h
@@ -916,6 +916,9 @@ public:
                      doublereal* phase_data);
 
 protected:
+    //! Cache for saved calculations within each Kinetics object.
+    ValueCache m_cache;
+
     // Update internal rate-of-progress variables #m_ropf and #m_ropr.
     virtual void updateROP() {
         throw NotImplementedError("Kinetics::updateROP");
diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp
index bdfaa0db8..40cab562a 100644
--- a/src/kinetics/InterfaceKinetics.cpp
+++ b/src/kinetics/InterfaceKinetics.cpp
@@ -184,6 +184,7 @@ void InterfaceKinetics::_update_rates_T()
 
         //  Calculate the forward rate constant by calling m_rates and store it in m_rfn[]
         m_rates.update(T, m_logtemp, DATA_PTR(m_rfn));
+        applyStickingCorrection(&m_rfn[0]);
 
         //  If we need to do conversions between exchange current density formulation and regular formulation
         //  (either way) do it here.
@@ -938,23 +939,34 @@ void InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
             }
         }
 
+        double surface_order = 0.0;
         // Adjust the A-factor
         for (Composition::const_iterator iter = r.reactants.begin();
              iter != r.reactants.end();
              ++iter) {
             size_t iPhase = speciesPhaseIndex(kineticsSpeciesIndex(iter->first));
             const ThermoPhase& p = thermo(iPhase);
+            const ThermoPhase& surf = thermo(surfacePhaseIndex());
             size_t k = p.speciesIndex(iter->first);
             if (iter->first == sticking_species) {
                 A_rate *= sqrt(GasConstant/(2*Pi*p.molecularWeight(k)));
             } else {
                 // Non-sticking species. Convert from coverages used in the
                 // sticking probability expression to the concentration units
-                // used in the mass action rate expression
+                // used in the mass action rate expression. For surface phases,
+                // the dependence on the site density is incorporated when the
+                // rate constant is evaluated, since we don't assume that the
+                // site density is known at this time.
                 double order = getValue(r.orders, iter->first, iter->second);
-                A_rate /= pow(p.standardConcentration(k), order);
+                if (&p == &surf) {
+                    A_rate *= pow(p.size(k), order);
+                    surface_order += order;
+                } else {
+                    A_rate *= pow(p.standardConcentration(k), -order);
+                }
             }
         }
+        m_sticking_orders.push_back(make_pair(m_ii, surface_order));
     }
     SurfaceArrhenius rate(A_rate, b_rate, r.rate.activationEnergy_R());
 
@@ -1314,6 +1326,30 @@ void InterfaceKinetics::determineFwdOrdersBV(ElectrochemicalReaction& r, std::ve
     }
 }
 
+void InterfaceKinetics::applyStickingCorrection(double* kf)
+{
+    if (m_sticking_orders.empty()) {
+        return;
+    }
+
+    static const int cacheId = m_cache.getId();
+    CachedArray cached = m_cache.getArray(cacheId);
+    vector_fp& factors = cached.value;
+
+    SurfPhase& surf = dynamic_cast<SurfPhase&>(thermo(reactionPhaseIndex()));
+    double n0 = surf.siteDensity();
+    if (!cached.validate(n0)) {
+        factors.resize(m_sticking_orders.size());
+        for (size_t n = 0; n < m_sticking_orders.size(); n++) {
+            factors[n] = pow(n0, -m_sticking_orders[n].second);
+        }
+    }
+
+    for (size_t n = 0; n < m_sticking_orders.size(); n++) {
+        kf[m_sticking_orders[n].first] *= factors[n];
+    }
+}
+
 
 void EdgeKinetics::finalize()
 {