From 259ceaa9a65261a3f56ccb2344c8f15a32727cd6 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Mon, 6 May 2013 20:05:53 +0000
Subject: [PATCH] Added assertions to check that some variables remain finite

---
 include/cantera/base/ctexceptions.h |  4 ++++
 src/kinetics/GasKinetics.cpp        | 11 +++++++++++
 src/zeroD/Reactor.cpp               | 10 ++++++++++
 3 files changed, 25 insertions(+)

diff --git a/include/cantera/base/ctexceptions.h b/include/cantera/base/ctexceptions.h
index 29093127d..dc5afcf81 100644
--- a/include/cantera/base/ctexceptions.h
+++ b/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
 
diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp
index 2a2483a43..c33de0546 100644
--- a/src/kinetics/GasKinetics.cpp
+++ b/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;
 }
 
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index 839432e8e..d5f32c349 100644
--- a/src/zeroD/Reactor.cpp
+++ b/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();