diff --git a/include/cantera/base/ctexceptions.h b/include/cantera/base/ctexceptions.h index 5cdb3eb3e..ede39719a 100644 --- a/include/cantera/base/ctexceptions.h +++ b/include/cantera/base/ctexceptions.h @@ -55,6 +55,37 @@ namespace Cantera * preprocessor symbol is defined, e.g. with the compiler option -DNDEBUG. */ +//! Enum containing Cantera's behavior for situations where overflow or underflow of real variables +//! may occur. +/*! + * Note this frequently occurs when taking exponentials of delta gibbs energies of reactions + * or when taking the exponentials of logs of activity coefficients. + */ +enum CT_RealNumber_Range_Behavior { + + //! For this specification of range behavior, nothing is done. This is the fastest + //! behavior when all calculations are believed to be ranged well. For situations + //! where there are range errors, NaN's or INF's will be introduced. + DONOTHING_CTRB = -1, + + //! For this specification of range behavior, the overflow or underflow calculation is changed. + //! Cantera will proceed by bounding the real number to maintain its viability, silently + //! changing the actual answer. + CHANGE_OVERFLOW_CTRB, + + //! When an overflow or underflow occurs, Cantera will throw an error + THROWON_OVERFLOW_CTRB, + + //! Cantera will use the fenv check capability introduced in C99 to check for + //! overflow and underflow conditions at crucial points. + //! It will throw an error if these conditions occur. + FENV_CHECK_CTRB, + + //! Cantera will throw an error in debug mode but will not in production mode. + //! (default) + THROWON_OVERFLOW_DEBUGMODEONLY_CTRB +}; + //! Base class for exceptions thrown by Cantera classes. /*! @@ -190,6 +221,16 @@ public: } }; +//! Quick check on whether there has been an underflow or overflow condition in the floating point unit +/*! + * @return Returns true if there has been such a condition and it has not been cleared. returns false + * if there hasn't been an overflow, underflow or invalid condition. + */ +extern bool check_FENV_OverUnder_Flow(); + +//! Clear all the flags for floating-point exceptions +extern void clear_FENV(); + //! Provides a line number #define XSTR_TRACE_LINE(s) STR_TRACE_LINE(s) diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index 7ee1fc392..c6016a343 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -7,6 +7,7 @@ #ifndef CT_PHASE_H #define CT_PHASE_H +#include "cantera/Cantera.h" #include "cantera/base/vec_functions.h" #include "cantera/base/ctml.h" #include "cantera/thermo/Elements.h" @@ -785,6 +786,14 @@ private: //! Entropy at 298.15 K and 1 bar of stable state pure elements (J kmol-1) vector_fp m_entropy298; +public: + //! Overflow behavior of real number calculations involving this thermo object + /*! + * The default is THROWON_OVERFLOW_CTRB + * Which throws an error in debug mode, but silently changes the answer in non-debug mode + */ + enum CT_RealNumber_Range_Behavior realNumberRangeBehavior_; + }; } diff --git a/src/base/ctexceptions.cpp b/src/base/ctexceptions.cpp index 2cafbb6e6..d77352195 100644 --- a/src/base/ctexceptions.cpp +++ b/src/base/ctexceptions.cpp @@ -5,6 +5,7 @@ #include "cantera/base/global.h" #include "cantera/base/stringUtils.h" +#include #include #include @@ -75,5 +76,19 @@ std::string IndexError::getMessage() const " outside valid range of 0 to " << (mmax_) << "."; return ss.str(); } +//============================================================================================================ +bool check_FENV_OverUnder_Flow() { + fexcept_t ff; + fegetexceptflag(&ff, FE_OVERFLOW || FE_UNDERFLOW || FE_INVALID); + if (ff) { + return true; + } + return false; +}; +//============================================================================================================ +void clear_FENV() { + feclearexcept(FE_ALL_EXCEPT); +} +//============================================================================================================ } // namespace Cantera diff --git a/src/thermo/GibbsExcessVPSSTP.cpp b/src/thermo/GibbsExcessVPSSTP.cpp index 063037563..eb3dcaaed 100644 --- a/src/thermo/GibbsExcessVPSSTP.cpp +++ b/src/thermo/GibbsExcessVPSSTP.cpp @@ -1,9 +1,8 @@ /** * @file GibbsExcessVPSSTP.cpp * Definitions for intermediate ThermoPhase object for phases which - * employ excess gibbs free energy formulations - * (see \ref thermoprops - * and class \link Cantera::GibbsExcessVPSSTP GibbsExcessVPSSTP\endlink). + * employ excess Gibbs free energy formulations + * (see \ref thermoprops and class \link Cantera::GibbsExcessVPSSTP GibbsExcessVPSSTP\endlink). * * Header file for a derived class of ThermoPhase that handles * variable pressure standard state methods for calculating @@ -25,7 +24,7 @@ using namespace std; namespace Cantera { - +//========================================================================================================================= GibbsExcessVPSSTP::GibbsExcessVPSSTP() : VPStandardStateTP(), moleFractions_(0), @@ -38,7 +37,7 @@ GibbsExcessVPSSTP::GibbsExcessVPSSTP() : m_pp(0) { } - +//========================================================================================================================= GibbsExcessVPSSTP::GibbsExcessVPSSTP(const GibbsExcessVPSSTP& b) : VPStandardStateTP(), moleFractions_(0), @@ -52,7 +51,7 @@ GibbsExcessVPSSTP::GibbsExcessVPSSTP(const GibbsExcessVPSSTP& b) : { GibbsExcessVPSSTP::operator=(b); } - +//========================================================================================================================= GibbsExcessVPSSTP& GibbsExcessVPSSTP::operator=(const GibbsExcessVPSSTP& b) { if (&b == this) { @@ -72,52 +71,52 @@ GibbsExcessVPSSTP& GibbsExcessVPSSTP::operator=(const GibbsExcessVPSSTP& b) return *this; } - +//========================================================================================================================= ThermoPhase* GibbsExcessVPSSTP::duplMyselfAsThermoPhase() const { return new GibbsExcessVPSSTP(*this); } - +//========================================================================================================================= void GibbsExcessVPSSTP::setMassFractions(const doublereal* const y) { Phase::setMassFractions(y); getMoleFractions(DATA_PTR(moleFractions_)); } - +//========================================================================================================================= void GibbsExcessVPSSTP::setMassFractions_NoNorm(const doublereal* const y) { Phase::setMassFractions_NoNorm(y); getMoleFractions(DATA_PTR(moleFractions_)); } - +//========================================================================================================================= void GibbsExcessVPSSTP::setMoleFractions(const doublereal* const x) { Phase::setMoleFractions(x); getMoleFractions(DATA_PTR(moleFractions_)); } - +//========================================================================================================================= void GibbsExcessVPSSTP::setMoleFractions_NoNorm(const doublereal* const x) { Phase::setMoleFractions_NoNorm(x); getMoleFractions(DATA_PTR(moleFractions_)); } - +//========================================================================================================================= void GibbsExcessVPSSTP::setConcentrations(const doublereal* const c) { Phase::setConcentrations(c); getMoleFractions(DATA_PTR(moleFractions_)); } - +//========================================================================================================================= /* * ------------ Mechanical Properties ------------------------------ */ - +//========================================================================================================================= void GibbsExcessVPSSTP::setPressure(doublereal p) { setState_TP(temperature(), p); } - +//========================================================================================================================= void GibbsExcessVPSSTP::calcDensity() { vector_fp vbar = getPartialMolarVolumesVector(); @@ -128,7 +127,7 @@ void GibbsExcessVPSSTP::calcDensity() doublereal dd = meanMolecularWeight() / vtotal; Phase::setDensity(dd); } - +//========================================================================================================================= void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p) { Phase::setTemperature(t); @@ -147,26 +146,25 @@ void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p) */ calcDensity(); } - +//========================================================================================================================= /* * - Activities, Standard States, Activity Concentrations ----------- */ - void GibbsExcessVPSSTP::getActivityConcentrations(doublereal* c) const { getActivities(c); } - +//========================================================================================================================= doublereal GibbsExcessVPSSTP::standardConcentration(size_t k) const { return 1.0; } - +//========================================================================================================================= doublereal GibbsExcessVPSSTP::logStandardConc(size_t k) const { return 0.0; } - +//========================================================================================================================= void GibbsExcessVPSSTP::getActivities(doublereal* ac) const { getActivityCoefficients(ac); @@ -175,24 +173,49 @@ void GibbsExcessVPSSTP::getActivities(doublereal* ac) const ac[k] *= moleFractions_[k]; } } - +//========================================================================================================================= void GibbsExcessVPSSTP::getActivityCoefficients(doublereal* const ac) const { - getLnActivityCoefficients(ac); - - // Protect against roundoff when taking exponentials - for (size_t k = 0; k < m_kk; k++) { - if (ac[k] > 700.) { - ac[k] = exp(700.0); - } else if (ac[k] < -700.) { - ac[k] = exp(-700.0); - } else { + // + // Protect against or inform about roundoff when taking exponentials + // + if ((DEBUG_MODE_ENABLED && realNumberRangeBehavior_ == THROWON_OVERFLOW_DEBUGMODEONLY_CTRB) || + (realNumberRangeBehavior_ == THROWON_OVERFLOW_CTRB)) { + for (size_t k = 0; k < m_kk; k++) { + if (ac[k] > 700.) { + throw CanteraError("GibbsExcessVPSSTP::getActivityCoefficients()", + "activity coefficient for " + int2str(k) + " is overflowing: ln(ac) = " + fp2str(ac[k])); + } else if (ac[k] < -700.) { + throw CanteraError("GibbsExcessVPSSTP::getActivityCoefficients()", + "activity coefficient for " + int2str(k) + " is underflowing: ln(ac) = " + fp2str(ac[k])); + } else { + ac[k] = exp(ac[k]); + } + } + } else if (realNumberRangeBehavior_ == CHANGE_OVERFLOW_CTRB) { + for (size_t k = 0; k < m_kk; k++) { + if (ac[k] > 700.) { + ac[k] = exp(700.0); + } else if (ac[k] < -700.) { + ac[k] = exp(-700.0); + } else { + ac[k] = exp(ac[k]); + } + } + } else { + for (size_t k = 0; k < m_kk; k++) { ac[k] = exp(ac[k]); } + if (realNumberRangeBehavior_ == FENV_CHECK_CTRB) { + if (check_FENV_OverUnder_Flow()) { + throw CanteraError("GibbsExcessVPSSTP::getActivityCoefficients()", + "activity coefficient is over/underflowing"); + } + } } } - +//========================================================================================================================= void GibbsExcessVPSSTP::getElectrochemPotentials(doublereal* mu) const { getChemPotentials(mu); @@ -201,11 +224,11 @@ void GibbsExcessVPSSTP::getElectrochemPotentials(doublereal* mu) const mu[k] += ve*charge(k); } } - +//========================================================================================================================= /* * ------------ Partial Molar Properties of the Solution ------------ */ - +//========================================================================================================================= void GibbsExcessVPSSTP::getPartialMolarVolumes(doublereal* vbar) const { /* @@ -213,12 +236,12 @@ void GibbsExcessVPSSTP::getPartialMolarVolumes(doublereal* vbar) const */ getStandardVolumes(vbar); } - +//========================================================================================================================= const vector_fp& GibbsExcessVPSSTP::getPartialMolarVolumesVector() const { return getStandardVolumes(); } - +//========================================================================================================================= double GibbsExcessVPSSTP::checkMFSum(const doublereal* const x) const { doublereal norm = accumulate(x, x + m_kk, 0.0); @@ -228,9 +251,13 @@ double GibbsExcessVPSSTP::checkMFSum(const doublereal* const x) const } return norm; } - +//========================================================================================================================= void GibbsExcessVPSSTP::getUnitsStandardConc(double* uA, int k, int sizeUA) const { + // + // We assume here that the units of the standard concentration is unitless. In other words activities are + // used unchanged in kinetics expressions. This may be changed in implementations of child classes. + // for (int i = 0; i < sizeUA; i++) { if (i == 0) { uA[0] = 0.0; @@ -252,14 +279,14 @@ void GibbsExcessVPSSTP::getUnitsStandardConc(double* uA, int k, int sizeUA) cons } } } - +//========================================================================================================================= void GibbsExcessVPSSTP::initThermo() { initLengths(); VPStandardStateTP::initThermo(); getMoleFractions(DATA_PTR(moleFractions_)); } - +//========================================================================================================================= void GibbsExcessVPSSTP::initLengths() { m_kk = nSpecies(); @@ -272,5 +299,5 @@ void GibbsExcessVPSSTP::initLengths() dlnActCoeffdlnN_.resize(m_kk, m_kk); m_pp.resize(m_kk); } - -} +//========================================================================================================================= +} // end of namespace Cantera diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp index ed96236bc..d4d45c951 100644 --- a/src/thermo/Phase.cpp +++ b/src/thermo/Phase.cpp @@ -27,7 +27,8 @@ Phase::Phase() : m_stateNum(-1), m_elementsFrozen(false), m_mm(0), - m_elem_type(0) + m_elem_type(0), + realNumberRangeBehavior_(THROWON_OVERFLOW_DEBUGMODEONLY_CTRB) { } @@ -43,7 +44,8 @@ Phase::Phase(const Phase& right) : m_stateNum(-1), m_elementsFrozen(false), m_mm(0), - m_elem_type(0) + m_elem_type(0), + realNumberRangeBehavior_(THROWON_OVERFLOW_DEBUGMODEONLY_CTRB) { // Use the assignment operator to do the actual copying operator=(right); @@ -106,6 +108,7 @@ Phase& Phase::operator=(const Phase& right) } m_id = right.m_id; m_name = right.m_name; + realNumberRangeBehavior_ = right.realNumberRangeBehavior_; return *this; }