From d76694119cea08720426f2bcf9585a6be98e4062 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 26 Oct 2007 15:34:49 +0000 Subject: [PATCH] Upgraded algorithm for setState_HP() and friends. Replaced the bare Newton's method with a root finder. This algorithm is tolerant of step jumps in Thermo parameters at the temperature boundaries. The algorithm also has an initial treatment of spinodals. Basically, if we are currently liquid, we want to stay liquid. If we are currently gas, we want to stay gas. --- Cantera/src/thermo/SpeciesThermoFactory.cpp | 6 +- Cantera/src/thermo/ThermoPhase.cpp | 535 +++++++++++++++++--- Cantera/src/thermo/ThermoPhase.h | 50 +- 3 files changed, 506 insertions(+), 85 deletions(-) diff --git a/Cantera/src/thermo/SpeciesThermoFactory.cpp b/Cantera/src/thermo/SpeciesThermoFactory.cpp index 835c6f2e3..10aba9a57 100755 --- a/Cantera/src/thermo/SpeciesThermoFactory.cpp +++ b/Cantera/src/thermo/SpeciesThermoFactory.cpp @@ -184,7 +184,7 @@ namespace Cantera { doublereal cplow = poly4(tmid, clow); doublereal cphigh = poly4(tmid, chigh); doublereal delta = cplow - cphigh; - if (fabs(delta/cplow) > 0.001) { + if (fabs(delta/(fabs(cplow)+1.0E-4)) > 0.001) { writelog("\n\n**** WARNING ****\nFor species "+name+ ", discontinuity in cp/R detected at Tmid = " +fp2str(tmid)+"\n"); @@ -198,7 +198,7 @@ namespace Cantera { doublereal hrtlow = enthalpy_RT(tmid, clow); doublereal hrthigh = enthalpy_RT(tmid, chigh); delta = hrtlow - hrthigh; - if (fabs(delta/hrtlow) > 0.001) { + if (fabs(delta/(fabs(hrtlow)+cplow*tmid)) > 0.001) { writelog("\n\n**** WARNING ****\nFor species "+name+ ", discontinuity in h/RT detected at Tmid = " +fp2str(tmid)+"\n"); @@ -212,7 +212,7 @@ namespace Cantera { doublereal srlow = entropy_R(tmid, clow); doublereal srhigh = entropy_R(tmid, chigh); delta = srlow - srhigh; - if (fabs(delta/srlow) > 0.001) { + if (fabs(delta/(fabs(srlow)+cplow)) > 0.001) { writelog("\n\n**** WARNING ****\nFor species "+name+ ", discontinuity in s/R detected at Tmid = " +fp2str(tmid)+"\n"); diff --git a/Cantera/src/thermo/ThermoPhase.cpp b/Cantera/src/thermo/ThermoPhase.cpp index e1c6bcd1b..0932675eb 100644 --- a/Cantera/src/thermo/ThermoPhase.cpp +++ b/Cantera/src/thermo/ThermoPhase.cpp @@ -22,6 +22,10 @@ #include "ThermoPhase.h" +#ifndef MAX +#define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif + using namespace std; namespace Cantera { @@ -191,85 +195,476 @@ namespace Cantera { setMassFractions(y); setPressure(p); } - void ThermoPhase::setState_HP(doublereal h, doublereal p, - doublereal tol) { - doublereal dt; - setPressure(p); + void ThermoPhase::setState_HP(doublereal Htarget, doublereal p, + doublereal dTtol) { + setState_HPorUV(Htarget, p, dTtol, false); + } - // Newton iteration - for (int n = 0; n < 500; n++) { - double h0 = enthalpy_mass(); - dt = (h - h0)/cp_mass(); - // limit step size to 100 K - if (dt > 100.0) dt = 100.0; - else if (dt < -100.0) dt = -100.0; - setState_TP(temperature() + dt, p); - if (fabs(dt) < tol) { - return; - } - } - throw CanteraError("setState_HP","No convergence. dt = " + fp2str(dt)); + void ThermoPhase::setState_UV(doublereal u, doublereal v, + doublereal dTtol) { + setState_HPorUV(u, v, dTtol, true); + } + + void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p, + doublereal dTtol, bool doUV) { + doublereal dt; + doublereal Hmax = 0.0, Hmin = 0.0;; + doublereal v = 0.0; + if (doUV) { + v = p; + setDensity(1.0/v); + } else { + setPressure(p); + } + double Tmax = maxTemp() + 0.1; + double Tmin = minTemp() - 0.1; + + // Make sure we are within the temperature bounds at the start + // of the iteration + double Tnew = temperature(); + if (Tnew > Tmax) { + Tnew = Tmax - 1.0; + if (doUV) { + setTemperature(Tnew); + } else { + setState_TP(Tnew, p); + } + } + if (Tnew < Tmin) { + Tnew = Tmin + 1.0; + if (doUV) { + setTemperature(Tnew); + } else { + setState_TP(Tnew, p); + } } - void ThermoPhase::setState_UV(doublereal u, doublereal v, - doublereal tol) { - doublereal dt; - setDensity(1.0/v); - for (int n = 0; n < 500; n++) { - dt = (u - intEnergy_mass())/cv_mass(); - if (dt > 100.0) dt = 100.0; - else if (dt < -100.0) dt = -100.0; - if (fabs(dt) < tol) { - setTemperature(temperature() + dt); - return; - } - setTemperature(temperature() + 0.5*dt); - } - throw CanteraError("setState_UV", - "no convergence. dt = " + fp2str(dt)+"\n" - +"tol = "+fp2str(tol)+"\n" - +"u = "+fp2str(u)+" v = "+fp2str(v)+"\n"); + double Hnew = 0.0; + double Cpnew = 0.0; + if (doUV) { + Hnew = intEnergy_mass(); + Cpnew = cv_mass(); + } else { + Hnew = enthalpy_mass(); + Cpnew = cp_mass(); + } + double Htop = Hnew; + double Ttop = Tnew; + double Hbot = Hnew; + double Tbot = Tnew; + double Told = Tnew; + double Hold = Hnew; + + bool ignoreBounds = false; + // Unstable phases are those for which + // cp < 0.0. These are possible for cases where + // we have passed the spinodal curve. + bool unstablePhase = false; + bool unstablePhaseNew = false; + + + // Newton iteration + for (int n = 0; n < 500; n++) { + Told = Tnew; + Hold = Hnew; + double cpd = Cpnew; + if (cpd < 0.0) { + 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 + if (dt > 100.0) dt = 100.0; + else if (dt < -100.0) dt = -100.0; + 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 + if (dt > 0.0) { + if (!unstablePhase) { + if (Htop > Htarget) { + if (Tnew > (0.75 * Ttop + 0.25 * Told)) { + dt = 0.75 * (Ttop - Told); + Tnew = Told + dt; + } + } + } else { + if (Hbot < Htarget) { + if (Tnew < (0.75 * Tbot + 0.25 * Told)) { + dt = 0.75 * (Tbot - Told); + Tnew = Told + dt; + } + } + } + } else { + if (!unstablePhase) { + if (Hbot < Htarget) { + if (Tnew < (0.75 * Tbot + 0.25 * Told)) { + dt = 0.75 * (Tbot - Told); + Tnew = Told + dt; + } + } + } else { + if (Htop > Htarget) { + if (Tnew > (0.75 * Ttop + 0.25 * Told)) { + dt = 0.75 * (Ttop - Told); + Tnew = Told + dt; + } + } + } + } + // Check Max and Min values + if (Tnew > Tmax) { + if (!ignoreBounds) { + if (doUV) { + setTemperature(Tmax); + Hmax = intEnergy_mass(); + } else { + setState_TP(Tmax, p); + Hmax = enthalpy_mass(); + } + if (Hmax >= Htarget) { + if (Htop < Htarget) { + Ttop = Tmax; + Htop = Hmax; + } + } else { + Tnew = Tmax + 1.0; + ignoreBounds = true; + } + } + } + if (Tnew < Tmin) { + if (!ignoreBounds) { + if (doUV) { + setTemperature(Tmin); + Hmin = intEnergy_mass(); + } else { + setState_TP(Tmin, p); + Hmin = enthalpy_mass(); + } + if (Hmin <= Htarget) { + if (Hbot > Htarget) { + Tbot = Tmin; + Hbot = Hmin; + } + } else { + Tnew = Tmin - 1.0; + ignoreBounds = true; + } + } + } + + // Try to keep phase within its region of stability + // -> Could do a lot better if I calculate the + // spinodal value of H. + for (int its = 0; its < 10; its++) { + Tnew = Told + dt; + if (doUV) { + setTemperature(Tnew); + Hnew = intEnergy_mass(); + Cpnew = cv_mass(); + } else { + setState_TP(Tnew, p); + Hnew = enthalpy_mass(); + Cpnew = cp_mass(); + } + if (Cpnew < 0.0) { + unstablePhaseNew = true; + } else { + break; + unstablePhaseNew = false; + } + if (unstablePhase == false) { + if (unstablePhaseNew == true) { + dt *= 0.25; + } + } + } + + if (Hnew == Htarget) { + return; + } else if (Hnew > Htarget) { + if ((Htop < Htarget) || (Hnew < Htop)) { + Htop = Hnew; + Ttop = Tnew; + } + } else if (Hnew < Htarget) { + if ((Hbot > Htarget) || (Hnew > Hbot)) { + Hbot = Hnew; + Tbot = Tnew; + } + } + // Convergence in H + double Herr = Htarget - Hnew; + double acpd = MAX(fabs(cpd), 1.0E-5); + double denom = MAX(fabs(Htarget), acpd * dTtol); + double HConvErr = fabs((Herr)/denom); + if (HConvErr < 0.00001 *dTtol) { + return; + } + if (fabs(dt) < dTtol) { + return; + } + } + throw CanteraError("setState_HPorUV","No convergence. dt = " + fp2str(dt)); + } + + void ThermoPhase::setState_SP(doublereal Starget, doublereal p, + doublereal dTtol) { + setState_SPorSV(Starget, p, dTtol, false); + } + + void ThermoPhase::setState_SV(doublereal Starget, doublereal v, + doublereal dTtol) { + setState_SPorSV(Starget, v, dTtol, true); + } + + void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p, + doublereal dTtol, bool doSV) { + doublereal v = 0.0; + doublereal dt; + if (doSV) { + v = p; + setDensity(1.0/v); + } else { + setPressure(p); + } + double Tmax = maxTemp() + 0.1; + double Tmin = minTemp() - 0.1; + + // Make sure we are within the temperature bounds at the start + // of the iteration + double Tnew = temperature(); + if (Tnew > Tmax) { + Tnew = Tmax - 1.0; + if (doSV) { + setTemperature(Tnew); + } else { + setState_TP(Tnew, p); + } + } + if (Tnew < Tmin) { + Tnew = Tmin + 1.0; + if (doSV) { + setTemperature(Tnew); + } else { + setState_TP(Tnew, p); + } } - void ThermoPhase::setState_SP(doublereal s, doublereal p, - doublereal tol) { - doublereal dt; - setPressure(p); - for (int n = 0; n < 500; n++) { - dt = (s - entropy_mass())*temperature()/cp_mass(); - if (dt > 100.0) dt = 100.0; - else if (dt < -100.0) dt = -100.0; - if (fabs(dt) < tol) { - setState_TP(temperature() + dt, p); - return; - } - setState_TP(temperature() + 0.5*dt, p); - } - throw CanteraError("setState_SP","no convergence. dt = " + fp2str(dt)); + double Snew = entropy_mass(); + double Cpnew = 0.0; + if (doSV) { + Cpnew = cv_mass(); + } else { + Cpnew = cp_mass(); } - void ThermoPhase::setState_SV(doublereal s, doublereal v, - doublereal tol) { - doublereal dt; - setDensity(1.0/v); - for (int n = 0; n < 500; n++) { - dt = (s - entropy_mass())*temperature()/cv_mass(); - if (dt > 100.0) dt = 100.0; - else if (dt < -100.0) dt = -100.0; - if (fabs(dt) < tol) { - setTemperature(temperature() + dt); - return; - } - setTemperature(temperature() + 0.5*dt); - } - throw CanteraError("setState_SV","no convergence. dt = " + fp2str(dt)); - } + double Stop = Snew; + double Ttop = Tnew; + double Sbot = Snew; + double Tbot = Tnew; + double Told = Tnew; + double Sold = Snew; - doublereal ThermoPhase::err(std::string msg) const { - throw CanteraError("ThermoPhase","Base class method " - +msg+" called. Equation of state type: "+int2str(eosType())); - return 0; + bool ignoreBounds = false; + // Unstable phases are those for which + // cp < 0.0. These are possible for cases where + // we have passed the spinodal curve. + bool unstablePhase = false; + bool unstablePhaseNew = false; + + + // Newton iteration + for (int n = 0; n < 500; n++) { + Told = Tnew; + Sold = Snew; + double cpd = Cpnew; + if (cpd < 0.0) { + 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; + else if (dt < -100.0) dt = -100.0; + Tnew = Told + dt; + // Limit the step size so that we are convergent + if (dt > 0.0) { + if (!unstablePhase) { + if (Stop > Starget) { + if (Tnew > Ttop) { + dt = 0.75 * (Ttop - Told); + Tnew = Told + dt; + } + } + } else { + if (Sbot < Starget) { + if (Tnew < Tbot) { + dt = 0.75 * (Tbot - Told); + Tnew = Told + dt; + } + } + } + } else { + if (!unstablePhase) { + if (Sbot < Starget) { + if (Tnew < Tbot) { + dt = 0.75 * (Tbot - Told); + Tnew = Told + dt; + } + } + } else { + if (Stop > Starget) { + if (Tnew > Ttop) { + dt = 0.75 * (Ttop - Told); + Tnew = Told + dt; + } + } + } + } + // Check Max and Min values + if (Tnew > Tmax) { + if (!ignoreBounds) { + if (doSV) { + setTemperature(Tmax); + } else { + setState_TP(Tmax, p); + } + double Smax = entropy_mass(); + if (Smax >= Starget) { + if (Stop < Starget) { + Ttop = Tmax; + Stop = Smax; + } + } else { + Tnew = Tmax + 1.0; + ignoreBounds = true; + } + } + } + if (Tnew < Tmin) { + if (!ignoreBounds) { + if (doSV) { + setTemperature(Tmin); + } else { + setState_TP(Tmin, p); + } + double Smin = enthalpy_mass(); + if (Smin <= Starget) { + if (Sbot > Starget) { + Sbot = Tmin; + Sbot = Smin; + } + } else { + Tnew = Tmin - 1.0; + ignoreBounds = true; + } + } + } + + // Try to keep phase within its region of stability + // -> Could do a lot better if I calculate the + // spinodal value of H. + for (int its = 0; its < 10; its++) { + Tnew = Told + dt; + if (doSV) { + setTemperature(Tnew); + Cpnew = cv_mass(); + } else { + setState_TP(Tnew, p); + Cpnew = cp_mass(); + } + Snew = entropy_mass(); + if (Cpnew < 0.0) { + unstablePhaseNew = true; + } else { + break; + unstablePhaseNew = false; + } + if (unstablePhase == false) { + if (unstablePhaseNew == true) { + dt *= 0.25; + } + } + } + + if (Snew == Starget) { + return; + } else if (Snew > Starget) { + if ((Stop < Starget) || (Snew < Stop)) { + Stop = Snew; + Ttop = Tnew; + } + } else if (Snew < Starget) { + if ((Sbot > Starget) || (Snew > Sbot)) { + Sbot = Snew; + Tbot = Tnew; + } + } + // Convergence in S + double Serr = Starget - Snew; + double acpd = MAX(fabs(cpd), 1.0E-5); + double denom = MAX(fabs(Starget), acpd * dTtol); + double SConvErr = fabs((Serr * Tnew)/denom); + if (SConvErr < 0.00001 *dTtol) { + return; + } + if (fabs(dt) < dTtol) { + return; + } } + throw CanteraError("setState_SPorSV","No convergence. dt = " + fp2str(dt)); + } + + doublereal ThermoPhase::err(std::string msg) const { + throw CanteraError("ThermoPhase","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + return 0; + } /* * Returns the units of the standard and general concentrations diff --git a/Cantera/src/thermo/ThermoPhase.h b/Cantera/src/thermo/ThermoPhase.h index 0c87a8b2d..296123393 100755 --- a/Cantera/src/thermo/ThermoPhase.h +++ b/Cantera/src/thermo/ThermoPhase.h @@ -1057,7 +1057,6 @@ namespace Cantera { */ void setState_PX(doublereal p, doublereal* x); - //! Set the internally storred pressure (Pa) and mass fractions. /*! * Note, the temperature is held constant during this operation. @@ -1074,11 +1073,10 @@ namespace Cantera { /*! * @param h Specific enthalpy (J/kg) * @param p Pressure (Pa) - * @param tol Optional parameter setting the tolerance of the - * calculation. + * @param tol Optional parameter setting the tolerance of the + * calculation. Defaults to 1.0E-4 */ - virtual void setState_HP(doublereal h, doublereal p, - doublereal tol = 1.e-4); + virtual void setState_HP(doublereal h, doublereal p, doublereal tol = 1.e-4); //! Set the specific internal energy (J/kg) and specific volume (m^3/kg). /*! @@ -1088,10 +1086,24 @@ namespace Cantera { * @param u specific internal energy (J/kg) * @param v specific volume (m^3/kg). * @param tol Optional parameter setting the tolerance of the - * calculation. + * calculation. Defaults to 1.0E-4 */ - virtual void setState_UV(doublereal u, doublereal v, - doublereal tol = 1.e-4); + virtual void setState_UV(doublereal u, doublereal v, doublereal tol = 1.e-4); + + private: + + //! Carry out work in HP and UV calculations. + /*! + * @param h Specific enthalpy or internal energy (J/kg) + * @param p Pressure (Pa) or specific volume (m^3/kg) + * @param tol Optional parameter setting the tolerance of the + * calculation. Defaults to 1.0E-4 + * @param doUV True if solving for UV, false for HP. + */ + void setState_HPorUV(doublereal h, doublereal p, + doublereal tol = 1.e-4, bool doUV = false); + + public: //! Set the specific entropy (J/kg/K) and pressure (Pa). /*! @@ -1101,10 +1113,9 @@ namespace Cantera { * @param s specific entropy (J/kg/K) * @param p specific pressure (Pa). * @param tol Optional parameter setting the tolerance of the - * calculation. + * calculation. Defaults to 1.0E-4 */ - virtual void setState_SP(doublereal s, doublereal p, - doublereal tol = 1.e-4); + virtual void setState_SP(doublereal s, doublereal p, doublereal tol = 1.e-4); //! Set the specific entropy (J/kg/K) and specific volume (m^3/kg). /*! @@ -1114,10 +1125,25 @@ namespace Cantera { * @param s specific entropy (J/kg/K) * @param v specific volume (m^3/kg). * @param tol Optional parameter setting the tolerance of the - * calculation. + * calculation. Defaults to 1.0E-4 */ virtual void setState_SV(doublereal s, doublereal v, doublereal tol = 1.e-4); + private: + + //! Carry out work in SP and SV calculations. + /*! + * @param s Specific entropy (J/kg) + * @param p Pressure (Pa) or specific volume (m^3/kg) + * @param tol Optional parameter setting the tolerance of the + * calculation. Defaults to 1.0E-4 + * @param doSV True if solving for SV, false for SP. + */ + void setState_SPorSV(doublereal s, doublereal p, + doublereal tol = 1.e-4, bool doSV = false); + + public: + //@} /**