research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Added in real number overflow exception control.

Browse source 
There are a series of activity coefficients for published battery models
 which actually overflows. Cantera was treating this as a nonerror and changing
the return results without comment. This is not correct. Made the behavior a
user controlled capability.
This commit is contained in:
Harry Moffat 2014-08-05 23:42:24 +00:00
parent cdec2e90c7
commit e129d02d93
5 changed files with 138 additions and 43 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

41
include/cantera/base/ctexceptions.h
View file

@ -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)

9
include/cantera/thermo/Phase.h
View file

@ -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_;
};
}

15
src/base/ctexceptions.cpp
View file

@ -5,6 +5,7 @@
#include "cantera/base/global.h"
#include "cantera/base/stringUtils.h"
#include <fenv.h>
#include <sstream>
#include <typeinfo>
@ -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

109
src/thermo/GibbsExcessVPSSTP.cpp
View file

@ -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

7
src/thermo/Phase.cpp
View file

@ -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;
}