From 6e138d0ae43f44f38eaee2b5c26f7e40b136265c Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 24 Aug 2015 12:27:07 -0400 Subject: [PATCH] [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 | 9 +++++++++ src/base/checkFinite.cpp | 18 ++++++++++++++++++ src/zeroD/IdealGasReactor.cpp | 10 ---------- src/zeroD/Reactor.cpp | 10 ---------- src/zeroD/ReactorNet.cpp | 2 ++ 5 files changed, 29 insertions(+), 20 deletions(-) diff --git a/include/cantera/base/utilities.h b/include/cantera/base/utilities.h index ff042173d..563444d89 100644 --- a/include/cantera/base/utilities.h +++ b/include/cantera/base/utilities.h @@ -521,6 +521,15 @@ void deepStdVectorPointerCopy(const std::vector &fromVec, std::vector &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' diff --git a/src/base/checkFinite.cpp b/src/base/checkFinite.cpp index 316e22757..134679eaf 100644 --- a/src/base/checkFinite.cpp +++ b/src/base/checkFinite.cpp @@ -13,6 +13,8 @@ #include #include +#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); + } + } +} + } diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 560a244c4..88852372b 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/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); } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 5e3f0a145..b47433ff9 100644 --- a/src/zeroD/Reactor.cpp +++ b/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); } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 8f429dae6..490fe9ce8 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/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]); }