[Reactor] Always check that input/output are finite

This causes problems to be caught early, where they can at least sometimes be
handled as recoverable intergrator errors, rather than letting NaN values
propagate through the solution.

include/cantera/base/utilities.h

@@ -521,6 +521,15 @@ void deepStdVectorPointerCopy(const std::vector<D*> &fromVec, std::vector<D*> &t
 //! Check to see that a number is finite (not NaN, +Inf or -Inf)
 void checkFinite(const double tmp);
 
+//! Check to see that all elements in an array are finite
+/*!
+ *  Throws an exception if any element is NaN, +Inf, or -Inf
+ *  @param name    Name to be used in the exception message if the check fails
+ *  @param values  Array of *N* values to be checked
+ *  @param N       Number of elements in *values*
+ */
+void checkFinite(const std::string& name, double* values, size_t N);
+
 //! Const accessor for a value in a std::map.
 /*!
  *  This is a const alternative to operator[]. Roughly equivalent to the 'at'

src/base/checkFinite.cpp

@@ -13,6 +13,8 @@
 
 #include <stdexcept>
 #include <cstdio>
+#include "cantera/base/stringUtils.h"
+#include "cantera/base/ctexceptions.h"
 
 // We expect that there will be special casing based on the computer
 // system here
@@ -54,4 +56,20 @@ void checkFinite(const double tmp)
     }
 }
 
+void checkFinite(const std::string& name, double* values, size_t N)
+{
+    for (size_t i = 0; i < N; i++) {
+        if (!finite(values[i])) {
+            std::string message = name + " contains non-finite elements:\n\n";
+            for (size_t j = 0; j < N; j++) {
+                if (!finite(values[j])) {
+                    message += name + "[" + int2str(j) + "] = " +
+                               fp2str(values[j]) + "\n";
+                }
+            }
+            throw CanteraError("checkFinite", message);
+        }
+    }
+}
+
 }

src/zeroD/IdealGasReactor.cpp

@@ -56,11 +56,6 @@ void IdealGasReactor::initialize(doublereal t0)
 
 void IdealGasReactor::updateState(doublereal* y)
 {
-    for (size_t i = 0; i < m_nv; i++) {
-        AssertFinite(y[i], "IdealGasReactor::updateState",
-                     "y[" + int2str(i) + "] is not finite");
-    }
-
     // The components of y are [0] the total mass, [1] the total volume,
     // [2] the temperature, [3...K+3] are the mass fractions of each species,
     // and [K+3...] are the coverages of surface species on each wall.
@@ -142,11 +137,6 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
         ydot[2] = 0;
     }
 
-    for (size_t i = 0; i < m_nv; i++) {
-        AssertFinite(ydot[i], "IdealGasReactor::evalEqs",
-                     "ydot[" + int2str(i) + "] is not finite");
-    }
-
     resetSensitivity(params);
 }
 

src/zeroD/Reactor.cpp

@@ -128,11 +128,6 @@ void Reactor::syncState()
 
 void Reactor::updateState(doublereal* y)
 {
-    for (size_t i = 0; i < m_nv; i++) {
-        AssertFinite(y[i], "Reactor::updateState",
-                     "y[" + int2str(i) + "] is not finite");
-    }
-
     // The components of y are [0] the total mass, [1] the total volume,
     // [2] the total internal energy, [3...K+3] are the mass fractions of each
     // species, and [K+3...] are the coverages of surface species on each wall.
@@ -267,11 +262,6 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
     }
 
     ydot[0] = dmdt;
-
-    for (size_t i = 0; i < m_nv; i++) {
-        AssertFinite(ydot[i], "Reactor::evalEqs",
-                     "ydot[" + int2str(i) + "] is not finite");
-    }
     resetSensitivity(params);
 }
 

src/zeroD/ReactorNet.cpp

@@ -156,6 +156,7 @@ void ReactorNet::eval(doublereal t, doublereal* y,
                                ydot + m_start[n], p + pstart);
         pstart += m_nparams[n];
     }
+    checkFinite("ydot", ydot, m_nv);
 }
 
 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
@@ -186,6 +187,7 @@ void ReactorNet::evalJacobian(doublereal t, doublereal* y,
 
 void ReactorNet::updateState(doublereal* y)
 {
+    checkFinite("y", y, m_nv);
     for (size_t n = 0; n < m_reactors.size(); n++) {
         m_reactors[n]->updateState(y + m_start[n]);
     }