[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

@@ -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) {

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: