diff --git a/include/cantera/base/mdp_allo.h b/include/cantera/base/mdp_allo.h index 39ac1d5fe..fc0fda045 100644 --- a/include/cantera/base/mdp_allo.h +++ b/include/cantera/base/mdp_allo.h @@ -687,46 +687,6 @@ extern void mdp_init_dbl_2(double** const v, const double value, */ extern void mdp_init_int_1(int* const v, const int value, const int len); - -/* - * Utility routines to check that a number is finite - */ - -//! Utility routine to check to see that a number is neither zero -//! nor indefinite. -/*! - * This check can be used before using the number in a denominator. - * - * @param tmp number to be checked - */ -extern void checkZeroFinite(const double tmp); - -//! Utility routine to check to see that a number is finite. -/*! - * @param tmp number to be checked - */ -extern void checkFinite(const double tmp); - -//! Utility routine to link checkFinte() to fortran program -/*! - * This routine is accessible from fortran, usually - * - * @param tmp Pointer to the number to check - * - * @todo link it into the usual way Cantera handles Fortran calls - */ -extern "C" void checkfinite_(double* tmp); - -//! utility routine to check that a double stays bounded -/*! - * This routine checks to see if a number stays bounded. The absolute - * value of the number is required to stay below the trigger. - * - * @param tmp Number to be checked - * @param trigger bounds on the number. Defaults to 1.0E20 - */ -extern void checkMagnitude(const double tmp, const double trigger = 1.0E20); - } /* end of mdp namespace */ /****************************************************************************/ #endif diff --git a/include/cantera/base/utilities.h b/include/cantera/base/utilities.h index f3930616c..d23ba6c6a 100644 --- a/include/cantera/base/utilities.h +++ b/include/cantera/base/utilities.h @@ -678,6 +678,9 @@ 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); } #endif diff --git a/src/base/checkFinite.cpp b/src/base/checkFinite.cpp index 5afb11dd1..5e6aa2c08 100644 --- a/src/base/checkFinite.cpp +++ b/src/base/checkFinite.cpp @@ -1,7 +1,6 @@ /** - * @file checkFinite.cpp - * Declarations for Multi Dimensional Pointer (mdp) routines that - * check for the presence of NaNs in the code. + * @file checkFinite.cpp Declarations for routines that check for the + * presence of NaNs in the code. */ /* * Copyright 2004 Sandia Corporation. Under the terms of Contract @@ -13,7 +12,6 @@ #include "cantera/base/ct_defs.h" #include -#include #include #include @@ -33,16 +31,11 @@ using namespace std; -namespace mdp -{ +namespace Cantera { -// Utility routine to check to see that a number is finite. -/* - * @param tmp number to be checked - */ -#ifdef _WIN32 void checkFinite(const double tmp) { +#ifdef _WIN32 if (_finite(tmp)) { if (_isnan(tmp)) { printf("checkFinite() ERROR: we have encountered a nan!\n"); @@ -51,14 +44,10 @@ void checkFinite(const double tmp) } else { printf("checkFinite() ERROR: we have encountered a neg inf!\n"); } - const std::string s = "checkFinite()"; - throw std::range_error(s); + throw std::range_error("checkFinite()"); } -} #else -void checkFinite(const double tmp) -{ - if (! finite(tmp)) { + if (!finite(tmp)) { if (isnan(tmp)) { printf("checkFinite() ERROR: we have encountered a nan!\n"); } else if (isinf(tmp) == 1) { @@ -66,90 +55,9 @@ void checkFinite(const double tmp) } else { printf("checkFinite() ERROR: we have encountered a neg inf!\n"); } - const std::string s = "checkFinite()"; - throw std::range_error(s); + throw std::range_error("checkFinite()"); } -} -#endif - - -// Utility routine to link checkFinte() to fortran program -/* - * This routine is accessible from fortran, usually - * - * @param tmp Pointer to the number to check - * - * @todo link it into the usual way Cantera handles Fortran calls - */ -extern "C" void checkfinite_(double* tmp) -{ - checkFinite(*tmp); -} - - -// Utility routine to check that a double stays bounded -/* - * This routine checks to see if a number stays bounded. The absolute - * value of the number is required to stay below the trigger. - * - * @param tmp Number to be checked - * @param trigger bounds on the number. Defaults to 1.0E20 - */ -void checkMagnitude(const double tmp, const double trigger) -{ - checkFinite(tmp); - if (fabs(tmp) >= trigger) { - char sbuf[64]; - sprintf(sbuf, "checkMagnitude() ERROR: Trigger %g exceeded: %g\n", trigger, - tmp); - throw std::range_error(sbuf); - } -} - -// Utility routine to check to see that a number is neither zero -// nor indefinite. -/* - * This check can be used before using the number in a denominator. - * - * @param tmp number to be checked - */ -#ifdef _WIN32 -void checkZeroFinite(const double tmp) -{ - if ((tmp == 0.0) || (! _finite(tmp))) { - if (tmp == 0.0) { - printf("checkZeroFinite() ERROR: we have encountered a zero!\n"); - } else if (_isnan(tmp)) { - printf("checkZeroFinite() ERROR: we have encountered a nan!\n"); - } else if (_fpclass(tmp) == _FPCLASS_PINF) { - printf("checkZeroFinite() ERROR: we have encountered a pos inf!\n"); - } else { - printf("checkZeroFinite() ERROR: we have encountered a neg inf!\n"); - } - char sbuf[64]; - sprintf(sbuf, "checkZeroFinite() ERROR: zero or indef exceeded: %g\n", - tmp); - throw std::range_error(sbuf); - } -} -#else -void checkZeroFinite(const double tmp) -{ - if ((tmp == 0.0) || (! finite(tmp))) { - if (tmp == 0.0) { - printf("checkZeroFinite() ERROR: we have encountered a zero!\n"); - } else if (isnan(tmp)) { - printf("checkZeroFinite() ERROR: we have encountered a nan!\n"); - } else if (isinf(tmp) == 1) { - printf("checkZeroFinite() ERROR: we have encountered a pos inf!\n"); - } else { - printf("checkZeroFinite() ERROR: we have encountered a neg inf!\n"); - } - char sbuf[64]; - sprintf(sbuf, "checkZeroFinite() ERROR: zero or indef exceeded: %g\n", - tmp); - throw std::range_error(sbuf); - } -} #endif } + +} diff --git a/src/numerics/NonlinearSolver.cpp b/src/numerics/NonlinearSolver.cpp index 1322aa674..5172015eb 100644 --- a/src/numerics/NonlinearSolver.cpp +++ b/src/numerics/NonlinearSolver.cpp @@ -574,17 +574,17 @@ doublereal NonlinearSolver::residErrorNorm(const doublereal* const resid, const for (size_t i = 0; i < neq_; i++) { #ifdef DEBUG_MODE - mdp::checkFinite(resid[i]); + checkFinite(resid[i]); #endif error = resid[i] / m_residWts[i]; #ifdef DEBUG_MODE - mdp::checkFinite(error); + checkFinite(error); #endif sum_norm += (error * error); } sum_norm = sqrt(sum_norm / neq_); #ifdef DEBUG_MODE - mdp::checkFinite(sum_norm); + checkFinite(sum_norm); #endif if (printLargest) { const int num_entries = printLargest; @@ -820,7 +820,7 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol]; } #ifdef DEBUG_MODE - mdp::checkFinite(m_rowWtScales[irow]); + checkFinite(m_rowWtScales[irow]); #endif jptr++; } @@ -844,7 +844,7 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, m_rowWtScales[irow] += vv * m_ewt[jcol]; } #ifdef DEBUG_MODE - mdp::checkFinite(m_rowWtScales[irow]); + checkFinite(m_rowWtScales[irow]); #endif } } @@ -1441,7 +1441,7 @@ doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac) deltaX_CP_[j] -= m_resid[i] * jac(i,j) * colFac * rowFac * m_ewt[j] * m_ewt[j] / (m_residWts[i] * m_residWts[i]); #ifdef DEBUG_MODE - mdp::checkFinite(deltaX_CP_[j]); + checkFinite(deltaX_CP_[j]); #endif } } @@ -3773,7 +3773,7 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, m_nJacEval++; #ifdef DEBUG_MODE for (int ii = 0; ii < neq_; ii++) { - mdp::checkFinite(f[ii]); + checkFinite(f[ii]); } #endif @@ -3844,7 +3844,7 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, throw CanteraError("NonlinearSolver::beuler_jac", "dy is equal to zero"); } for (int ii = 0; ii < neq_; ii++) { - mdp::checkFinite(m_wksp[ii]); + checkFinite(m_wksp[ii]); } #endif @@ -3929,7 +3929,7 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, throw CanteraError("NonlinearSolver::beuler_jac", "dy is equal to zero"); } for (int ii = 0; ii < neq_; ii++) { - mdp::checkFinite(m_wksp[ii]); + checkFinite(m_wksp[ii]); } #endif if (info != 1) { @@ -4094,7 +4094,7 @@ NonlinearSolver::computeResidWts() for (size_t i = 0; i < neq_; i++) { m_residWts[i] = userResidAtol_[i] + userResidRtol_ * m_rowWtScales[i] / rtol_; #ifdef DEBUG_MODE - mdp::checkFinite(m_residWts[i]); + checkFinite(m_residWts[i]); #endif } } else { @@ -4102,7 +4102,7 @@ NonlinearSolver::computeResidWts() for (size_t i = 0; i < neq_; i++) { m_residWts[i] = m_rowWtScales[i]; #ifdef DEBUG_MODE - mdp::checkFinite(m_residWts[i]); + checkFinite(m_residWts[i]); #endif sum += m_residWts[i]; } diff --git a/src/numerics/RootFind.cpp b/src/numerics/RootFind.cpp index 51bd5be51..981b87dc6 100644 --- a/src/numerics/RootFind.cpp +++ b/src/numerics/RootFind.cpp @@ -18,9 +18,7 @@ #endif #include "cantera/base/global.h" -#ifdef DEBUG_MODE -#include "cantera/base/mdp_allo.h" -#endif +#include "cantera/base/utilities.h" #include "cantera/base/stringUtils.h" /* Standard include files */ @@ -1183,11 +1181,11 @@ doublereal RootFind::func(doublereal x) { doublereal r; #ifdef DEBUG_MODE - mdp::checkFinite(x); + checkFinite(x); #endif m_residFunc->evalSS(0.0, &x, &r); #ifdef DEBUG_MODE - mdp::checkFinite(r); + checkFinite(r); #endif doublereal ff = r - m_funcTargetValue; if (x >= x_maxTried_) {