[Equil] Avoid floating point overflow from pressure term

Since computing the residual norm effectively requires p**2, the
correct limit to avoid an overflow is sqrt(DBL_MAX), or roughly
exp(300).

This change fixes crashes due to SIGFPE in cases where floating point
exceptions have been enabled, and also fixes a few cases where the
solver did not previously converge.
This commit is contained in:
Ray Speth 2013-09-30 01:53:54 +00:00
parent 2ac98393bd
commit f95bc8fb43

View file

@ -380,10 +380,10 @@ void IdealGasPhase::setToEquilState(const doublereal* mu_RT)
tmp = -grt[k] + mu_RT[k];
if (tmp < -600.) {
m_pp[k] = 0.0;
} else if (tmp > 500.0) {
tmp2 = tmp / 500.;
} else if (tmp > 300.0) {
tmp2 = tmp / 300.;
tmp2 *= tmp2;
m_pp[k] = m_p0 * exp(500.) * tmp2;
m_pp[k] = m_p0 * exp(300.) * tmp2;
} else {
m_pp[k] = m_p0 * exp(tmp);
}