From a6a10b854c552ed48ace40cbd1d8103e9ac8fc5c Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 7 Nov 2014 02:14:55 +0000 Subject: [PATCH] [Reactor] Add damping factor to setState Newton iteration This resolves convergence issues in somes cases where the iteration would oscillate around the true solution. --- src/zeroD/Reactor.cpp | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index d8b3223a7..bcd91852f 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -137,20 +137,32 @@ void Reactor::updateState(doublereal* y) m_thermo->setMassFractions_NoNorm(y+3); if (m_energy) { - // Use Newton's method to determine the mixture temperature. Tight - // tolerances are required both for Jacobian evaluation and for + // Use a damped Newton's method to determine the mixture temperature. + // Tight tolerances are required both for Jacobian evaluation and for // sensitivity analysis to work correctly. doublereal U = y[2]; doublereal T = temperature(); double dT = 100; + double dUprev = 1e10; + double dU = 1e10; int i = 0; + double damp = 1.0; while (abs(dT / T) > 10 * DBL_EPSILON) { + dUprev = dU; m_thermo->setState_TR(T, m_mass / m_vol); double dUdT = m_thermo->cv_mass() * m_mass; - dT = (m_thermo->intEnergy_mass() * m_mass - U) / dUdT; - dT = std::min(dT, 0.5 * T); + dU = m_thermo->intEnergy_mass() * m_mass - U; + dT = dU / dUdT; + // Reduce the damping coefficient if the magnitude of the error + // isn't decreasing + if (std::abs(dU) < std::abs(dUprev)) { + damp = 1.0; + } else { + damp *= 0.8; + } + dT = std::min(dT, 0.5 * T) * damp; T -= dT; i++; if (i > 100) {