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

Fixed an error in setState_UV() and similar routines that was occasionally

Browse source 
causing numerical errors in the calculation of Jacobians.
If the tolerances to these routines were set too low for calculation
of Jacobians, then the dt*=1.5 logic was causing the evaluation to
kick out of the routine with an incorrect deltaU value on the first
iteration. This meant that the numerically derived Jacobian entry
was off by a factor of 50%.
This commit is contained in:
Harry Moffat 2007-12-04 17:34:19 +00:00
parent 5060afcec3
commit 235958d810
1 changed files with 9 additions and 48 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

57
Cantera/src/thermo/ThermoPhase.cpp
View file

@ -62,12 +62,12 @@ namespace Cantera {
m_hasElementPotentials(false),
m_chargeNeutralityNecessary(false)
{
/*
* Call the assignment operator
*/
*this = operator=(right);
/*
* Call the assignment operator
*/
*this = operator=(right);
}
/*
* operator=()
*
@ -272,32 +272,14 @@ namespace Cantera {
unstablePhase = true;
}
dt = (Htarget - Hold)/cpd;
if (dt > 0.0) {
if (!unstablePhase) {
if (Htop < Htarget) {
dt *= 1.5;
}
} else {
if (Hbot > Htarget) {
dt *= 1.5;
}
}
} else {
if (!unstablePhase) {
if (Hbot > Htarget) {
dt *= 1.5;
}
} else {
if (Htop < Htarget) {
dt *= 1.5;
}
}
}
// limit step size to 200 K
// limit step size to 210 K
if (dt > 100.0) dt = 100.0;
else if (dt < -100.0) dt = -100.0;
// Calculate the new T
Tnew = Told + dt;
// Limit the step size so that we are convergent
// This is the step that makes it different from a
// Newton's algorithm
@ -506,27 +488,6 @@ namespace Cantera {
unstablePhase = true;
}
dt = (Starget - Sold)*Told/cpd;
if (dt > 0.0) {
if (!unstablePhase) {
if (Stop < Starget) {
dt *= 1.5;
}
} else {
if (Sbot > Starget) {
dt *= 1.5;
}
}
} else {
if (!unstablePhase) {
if (Sbot > Starget) {
dt *= 1.5;
}
} else {
if (Stop < Starget) {
dt *= 1.5;
}
}
}
// limit step size to 200 K
if (dt > 100.0) dt = 100.0;