Added assertions to check that some variables remain finite

include/cantera/base/ctexceptions.h

@@ -248,6 +248,10 @@ private:
 #endif
 
 
+//! Throw an exception if the specified exception is not a finite number.
+#ifndef AssertFinite
+#  define AssertFinite(expr, procedure, message) AssertThrowMsg(expr < BigNumber && expr > -BigNumber, procedure, message)
+#endif
 
 #endif
 

src/kinetics/GasKinetics.cpp

@@ -341,6 +341,8 @@ void GasKinetics::processFalloffReactions()
 
     for (size_t i = 0; i < m_nfall; i++) {
         pr[i] = concm_falloff_values[i] * m_rfn_low[i] / m_rfn_high[i];
+        AssertFinite(pr[i], "GasKinetics::processFalloffReactions",
+                     "pr[" + int2str(i) + "] is not finite.");
     }
 
     double* work = (falloff_work.empty()) ? 0 : &falloff_work[0];
@@ -399,6 +401,15 @@ void GasKinetics::updateROP()
         m_ropnet[j] = m_ropf[j] - m_ropr[j];
     }
 
+    for (size_t i = 0; i < m_rfn.size(); i++) {
+        AssertFinite(m_rfn[i], "GasKinetics::updateROP",
+                     "m_rfn[" + int2str(i) + "] is not finite.");
+        AssertFinite(m_ropf[i], "GasKinetics::updateROP",
+                     "m_ropf[" + int2str(i) + "] is not finite.");
+        AssertFinite(m_ropr[i], "GasKinetics::updateROP",
+                     "m_ropr[" + int2str(i) + "] is not finite.");
+    }
+
     m_ROP_ok = true;
 }
 

src/zeroD/Reactor.cpp

@@ -125,6 +125,11 @@ size_t Reactor::nSensParams()
 
 void Reactor::updateState(doublereal* y)
 {
+    for (size_t i = 0; i < m_nsp+2; i++) {
+        AssertFinite(y[i], "Reactor::updateState",
+                     "y[" + int2str(i) + "] is not finite");
+    }
+
     // The components of y are [0] the total internal energy,
     // [1] the total volume, and [2...K+2] the mass of each species.
     m_vol = y[1];
@@ -296,6 +301,11 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
         }
     }
 
+    for (size_t i = 0; i < m_nsp+2; i++) {
+        AssertFinite(ydot[i], "Reactor::evalEqs",
+                     "ydot[" + int2str(i) + "] is not finite");
+    }
+
     // reset sensitivity parameters
     if (params) {
         size_t npar = m_pnum.size();