Update importKinetics to identify unspecified electrochemical reactions

Add test coverage for beta default value for electrochem reactions

include/cantera/kinetics/importKinetics.h

@@ -110,6 +110,26 @@ bool importKinetics(const XML_Node& phase, std::vector<ThermoPhase*> th,
  */
 bool buildSolutionFromXML(XML_Node& root, const std::string& id,
                           const std::string& nm, ThermoPhase* th, Kinetics* kin);
+
+//! Check to ensure that all electrochemical reactions are specified correctly
+/*!
+ *  This function ensures the user has correctly specified all electrochemical
+ *  reactions. The routine counts the amount of charge (i.e. number of electron
+ *  elements specified for each species in each phase) for both reactants and
+ *  products. If net charge transfer phases during a reaction, the reaction is
+ *  electrochemical. If not already specified as such, the function defines the
+ *  reaction as electrochemical, corrects the reaction attributes, and sets
+ *  beta = 0.5.
+ *
+ * @param p     This is an XML node containing a description of the owning
+ *              phase for the kinetics object.
+ * @param kin   This is a pointer to a kinetics manager class.
+ * @param r     This is the reaction node that is being evaluated
+ * @return      The function always returns true.
+*/
+bool checkElectrochemReaction(const XML_Node& p, Kinetics& kin, const XML_Node& r);
+
+
 }
 
 #endif

src/kinetics/Reaction.cpp

@@ -255,7 +255,7 @@ InterfaceReaction::InterfaceReaction(const Composition& reactants_,
 
 ElectrochemicalReaction::ElectrochemicalReaction()
     : film_resistivity(0.0)
-    , beta(0.0)
+    , beta(0.5)
     , exchange_current_density_formulation(false)
 {
 }
@@ -265,12 +265,11 @@ ElectrochemicalReaction::ElectrochemicalReaction(const Composition& reactants_,
                                                  const Arrhenius& rate_)
     : InterfaceReaction(reactants_, products_, rate_)
     , film_resistivity(0.0)
-    , beta(0.0)
+    , beta(0.5)
     , exchange_current_density_formulation(false)
 {
 }
 
-
 Arrhenius readArrhenius(const XML_Node& arrhenius_node)
 {
     return Arrhenius(getFloat(arrhenius_node, "A", "toSI"),

src/kinetics/importKinetics.cpp

@@ -79,6 +79,7 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin,
         // if no 'include' directive, then include all reactions
         if (incl.empty()) {
             for (size_t i = 0; i < allrxns.size(); i++) {
+                checkElectrochemReaction(p,kin,*allrxns[i]);
                 kin.addReaction(newReaction(*allrxns[i]));
                 ++itot;
             }
@@ -110,6 +111,7 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin,
                         // do a lexical min max and operation. This sometimes
                         // has surprising results.
                         if ((rxid >= imin) && (rxid <= imax)) {
+                            checkElectrochemReaction(p,kin,*r);
                             kin.addReaction(newReaction(*r));
                             ++itot;
                         }
@@ -146,7 +148,7 @@ bool importKinetics(const XML_Node& phase, std::vector<ThermoPhase*> th,
         }
     }
 
-    // if other phases are involved in the reaction mechanism, they must be
+    // If other phases are involved in the reaction mechanism, they must be
     // listed in a 'phaseArray' child element. Homogeneous mechanisms do not
     // need to include a phaseArray element.
     vector<string> phase_ids;
@@ -212,4 +214,87 @@ bool buildSolutionFromXML(XML_Node& root, const std::string& id,
     return true;
 }
 
+bool checkElectrochemReaction(const XML_Node& p, Kinetics& kin, const XML_Node& r)
+{
+    // If other phases are involved in the reaction mechanism, they must be
+    // listed in a 'phaseArray' child element. Homogeneous mechanisms do not
+    // need to include a phaseArray element.
+    vector<string> phase_ids;
+    if (p.hasChild("phaseArray")) {
+        const XML_Node& pa = p.child("phaseArray");
+        getStringArray(pa, phase_ids);
+    }
+    phase_ids.push_back(p["id"]);
+
+    // Get reaction product and reactant information
+    Composition reactants = parseCompString(r.child("reactants").value());
+    Composition products = parseCompString(r.child("products").value());
+
+
+    // If the reaction has undeclared species don't perform electrochemical check
+    for (const auto& sp : reactants) {
+        if (kin.kineticsSpeciesIndex(sp.first) == npos) {
+            return true;
+        }
+    }
+
+    for (const auto& sp : products) {
+        if (kin.kineticsSpeciesIndex(sp.first) == npos) {
+            return true;
+        }
+    }
+
+    // Initialize the electron counter for each phase
+    std::vector<double> e_counter(phase_ids.size(), 0.0);
+
+    // Find the amount of electrons in the products for each phase
+    for (const auto& sp : products) {
+        const ThermoPhase& ph = kin.speciesPhase(sp.first);
+        size_t k = ph.speciesIndex(sp.first);
+        double stoich = sp.second;
+        for (size_t m = 0; m < phase_ids.size(); m++) {
+            if (phase_ids[m] == ph.id()) {
+                e_counter[m] += stoich * ph.charge(k);
+                break;
+            }
+        }
+    }
+
+    // Subtract the amount of electrons in the reactants for each phase
+    for (const auto& sp : reactants) {
+        const ThermoPhase& ph = kin.speciesPhase(sp.first);
+        size_t k = ph.speciesIndex(sp.first);
+        double stoich = sp.second;
+        for (size_t m = 0; m < phase_ids.size(); m++) {
+            if (phase_ids[m] == ph.id()) {
+                e_counter[m] -= stoich * ph.charge(k);
+                break;
+            }
+        }
+    }
+
+    // If the electrons change phases then the reaction is electrochemical
+    bool echemical = false;
+    for(size_t m = 0; m < phase_ids.size(); m++) {
+        if (fabs(e_counter[m]) > 1e-4) {
+            echemical = true;
+            break;
+        }
+    }
+
+    // If the reaction is electrochemical, ensure the reaction is identified as
+    // electrochemical. If not already specified beta is assumed to be 0.5
+    std::string type = ba::to_lower_copy(r["type"]);
+    if (!r.child("rateCoeff").hasChild("electrochem")) {
+        if ((type != "butlervolmer_noactivitycoeffs" &&
+             type != "butlervolmer" &&
+             type != "surfaceaffinity") &&
+             echemical) {
+            XML_Node& f = r.child("rateCoeff").addChild("electrochem","");
+            f.addAttribute("beta",0.5);
+        }
+    }
+    return true;
+}
+
 }

test/data/sofc-test.xml

@@ -412,7 +412,6 @@
     <reaction reversible="yes" type="edge" id="edge-f2">
       <equation>H(m) + O''(ox) [=] (m) + electron + OH'(ox)</equation>
       <rateCoeff>
-        <electrochem beta="0.5"/>
         <Arrhenius>
            <A>5.000000E+10</A>
            <b>0.0</b>