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: + //@} /**