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

[Reactor] Add damping factor to setState Newton iteration

Browse source 
This resolves convergence issues in somes cases where the iteration would
oscillate around the true solution.
This commit is contained in:
Ray Speth 2014-11-07 02:14:55 +00:00
parent 74a522bdda
commit a6a10b854c
1 changed files with 16 additions and 4 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

20
src/zeroD/Reactor.cpp
View file

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