From ac5337130a8d21728026f04155c0c85ad2ea0c5e Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Fri, 20 Jan 2017 17:21:42 -0500
Subject: [PATCH] [Kinetics] Validate balance of surface sites for interface
 reactions

The number of surface sites should be the same for the reactants and products.

Fixes #412
---
 src/kinetics/InterfaceKinetics.cpp    | 25 +++++++++++++++++++++++++
 test/kinetics/kineticsFromScratch.cpp | 10 ++++++++++
 2 files changed, 35 insertions(+)

diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp
index 8418feb17..1d5a617a6 100644
--- a/src/kinetics/InterfaceKinetics.cpp
+++ b/src/kinetics/InterfaceKinetics.cpp
@@ -604,6 +604,31 @@ void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
 
 bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
 {
+    if (!m_surf) {
+        init();
+    }
+
+    // Check that the number of surface sites is balanced
+    double reac_sites = 0.0;
+    double prod_sites = 0.0;
+    for (const auto& reactant : r_base->reactants) {
+        size_t k = m_surf->speciesIndex(reactant.first);
+        if (k != npos) {
+            reac_sites += reactant.second * m_surf->size(k);
+        }
+    }
+    for (const auto& product : r_base->products) {
+        size_t k = m_surf->speciesIndex(product.first);
+        if (k != npos) {
+            prod_sites += product.second * m_surf->size(k);
+        }
+    }
+    if (fabs(reac_sites - prod_sites) > 1e-5 * (reac_sites + prod_sites)) {
+        throw CanteraError("InterfaceKinetics::addReaction", "Number of surface"
+            " sites not balanced in reaction {}.\nReactant sites: {}\n"
+            "Product sites: {}", r_base->equation(), reac_sites, prod_sites);
+    }
+
     size_t i = nReactions();
     bool added = Kinetics::addReaction(r_base);
     if (!added) {
diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp
index 22ecaa8b2..74e7cbf85 100644
--- a/test/kinetics/kineticsFromScratch.cpp
+++ b/test/kinetics/kineticsFromScratch.cpp
@@ -379,6 +379,16 @@ TEST_F(InterfaceKineticsFromScratch, add_sticking_reaction)
     check_rates(0);
 }
 
+TEST_F(InterfaceKineticsFromScratch, unbalanced_sites)
+{
+    Composition reac = parseCompString("H(m):1 O(m):1");
+    Composition prod = parseCompString("OH(m):1");
+    Arrhenius rate(5e21, 0, 100.0e6 / GasConstant);
+
+    auto R = make_shared<InterfaceReaction>(reac, prod, rate);
+    ASSERT_THROW(kin.addReaction(R), CanteraError);
+}
+
 class KineticsAddSpecies : public testing::Test
 {
 public: