From d72c2e7132ea0069b7e6c421d27224ccf9c2fe0f Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Sat, 15 Nov 2014 00:46:57 +0000
Subject: [PATCH] [Kinetics] Add reaction balance check to Kinetics

This applies to reactions added using Reaction objects. Previously, the
check would only be performed for reactions added from an XML tree.
---
 include/cantera/kinetics/Kinetics.h |  5 +++
 src/kinetics/Kinetics.cpp           | 48 +++++++++++++++++++++++++++++
 2 files changed, 53 insertions(+)

diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h
index 93e34f021..b2076f4d6 100644
--- a/include/cantera/kinetics/Kinetics.h
+++ b/include/cantera/kinetics/Kinetics.h
@@ -915,6 +915,11 @@ protected:
     double checkDuplicateStoich(std::map<int, double>& r1,
                                 std::map<int, double>& r2) const;
 
+    //! Check that the specified reaction is balanced (same number of atoms for
+    //! each element in the reactants and products). Raises an exception if the
+    //! reaction is not balanced.
+    void checkReactionBalance(const Reaction& R);
+
     //! @name Stoichiometry management
     /*!
      *  These objects and functions handle turning reaction extents into species
diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp
index df9baac8c..b56b394ee 100644
--- a/src/kinetics/Kinetics.cpp
+++ b/src/kinetics/Kinetics.cpp
@@ -288,6 +288,52 @@ double Kinetics::checkDuplicateStoich(std::map<int, double>& r1,
     return ratio;
 }
 
+void Kinetics::checkReactionBalance(const Reaction& R)
+{
+    Composition balr, balp;
+    // iterate over the products
+    for (Composition::const_iterator iter = R.products.begin();
+         iter != R.products.end();
+         ++iter) {
+        const ThermoPhase& ph = speciesPhase(iter->first);
+        size_t k = ph.speciesIndex(iter->first);
+        double stoich = iter->second;
+        for (size_t m = 0; m < ph.nElements(); m++) {
+            balr[ph.elementName(m)] = 0.0; // so that balr contains all species
+            balp[ph.elementName(m)] += stoich*ph.nAtoms(k,m);
+        }
+    }
+    for (Composition::const_iterator iter = R.reactants.begin();
+         iter != R.reactants.end();
+         ++iter) {
+        const ThermoPhase& ph = speciesPhase(iter->first);
+        size_t k = ph.speciesIndex(iter->first);
+        double stoich = iter->second;
+        for (size_t m = 0; m < ph.nElements(); m++) {
+            balr[ph.elementName(m)] += stoich*ph.nAtoms(k,m);
+        }
+    }
+
+    string msg;
+    bool ok = true;
+    for (Composition::iterator iter = balr.begin();
+         iter != balr.end();
+         ++iter) {
+        const string& elem = iter->first;
+        double elemsum = balr[elem] + balp[elem];
+        double elemdiff = fabs(balp[elem] - balr[elem]);
+        if (elemsum > 0.0 && elemdiff/elemsum > 1e-5) {
+            ok = false;
+            msg += "  " + elem + "           " + fp2str(balr[elem]) +
+                   "           " + fp2str(balp[elem]) + "\n";
+        }
+    }
+    if (!ok) {
+        msg = "The following reaction is unbalanced: " + R.equation() + "\n" +
+              "  Element    Reactants    Products\n" + msg;
+        throw CanteraError("checkRxnElementBalance", msg);
+    }
+}
 
 void Kinetics::selectPhase(const doublereal* data, const thermo_t* phase,
                            doublereal* phase_data)
@@ -602,6 +648,8 @@ void Kinetics::addReaction(ReactionData& r) {
 
 void Kinetics::addReaction(shared_ptr<Reaction> r)
 {
+    checkReactionBalance(*r);
+
     size_t irxn = nReactions();
     std::vector<size_t> rk, pk;
     vector_fp rstoich, pstoich;