From f95bc8fb4381d50842d0c51c85e28af40805eb3d Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 30 Sep 2013 01:53:54 +0000 Subject: [PATCH] [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. --- src/thermo/IdealGasPhase.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/thermo/IdealGasPhase.cpp b/src/thermo/IdealGasPhase.cpp index 9d7519b28..5fbde3100 100644 --- a/src/thermo/IdealGasPhase.cpp +++ b/src/thermo/IdealGasPhase.cpp @@ -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); }