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

[Thermo] Make tolerances used in state setters more consistent

Browse source
This commit is contained in:
Ray Speth 2016-10-18 18:10:18 -04:00
parent 5b1a4a60bb
commit 478a62d2af
6 changed files with 71 additions and 79 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
include/cantera/thermo/PureFluidPhase.h
View file

@ -129,17 +129,17 @@ public:
* @{ * @{
*/ */
virtual void setState_HP(double h, double p, double tol = 1.e-8); virtual void setState_HP(double h, double p, double tol=1e-9);
virtual void setState_UV(double u, double v, double tol = 1.e-8); virtual void setState_UV(double u, double v, double tol=1e-9);
virtual void setState_SV(double s, double v, double tol = 1.e-8); virtual void setState_SV(double s, double v, double tol=1e-9);
virtual void setState_SP(double s, double p, double tol = 1.e-8); virtual void setState_SP(double s, double p, double tol=1e-9);
virtual void setState_ST(double s, double t, double tol = 1.e-8); virtual void setState_ST(double s, double t, double tol=1e-9);
virtual void setState_TV(double t, double v, double tol = 1.e-8); virtual void setState_TV(double t, double v, double tol=1e-9);
virtual void setState_PV(double p, double v, double tol = 1.e-8); virtual void setState_PV(double p, double v, double tol=1e-9);
virtual void setState_UP(double u, double p, double tol = 1.e-8); virtual void setState_UP(double u, double p, double tol=1e-9);
virtual void setState_VH(double v, double h, double tol = 1.e-8); virtual void setState_VH(double v, double h, double tol=1e-9);
virtual void setState_TH(double t, double h, double tol = 1.e-8); virtual void setState_TH(double t, double h, double tol=1e-9);
virtual void setState_SH(double s, double h, double tol = 1.e-8); virtual void setState_SH(double s, double h, double tol=1e-9);
//@} //@}
//! @name Critical State Properties //! @name Critical State Properties

12
include/cantera/thermo/SingleSpeciesTP.h
View file

@ -245,14 +245,10 @@ public:
//! Mole fractions are fixed, with x[0] = 1.0. //! Mole fractions are fixed, with x[0] = 1.0.
virtual void setMoleFractions(const doublereal* const x) {}; virtual void setMoleFractions(const doublereal* const x) {};
virtual void setState_HP(doublereal h, doublereal p, virtual void setState_HP(double h, double p, double tol=1e-9);
doublereal tol = 1.e-8); virtual void setState_UV(double u, double v, double tol=1e-9);
virtual void setState_UV(doublereal u, doublereal v, virtual void setState_SP(double s, double p, double tol=1e-9);
doublereal tol = 1.e-8); virtual void setState_SV(double s, double v, double tol=1e-9);
virtual void setState_SP(doublereal s, doublereal p,
doublereal tol = 1.e-8);
virtual void setState_SV(doublereal s, doublereal v,
doublereal tol = 1.e-8);
//@} //@}
virtual bool addSpecies(shared_ptr<Species> spec); virtual bool addSpecies(shared_ptr<Species> spec);

75
include/cantera/thermo/ThermoPhase.h
View file

@ -948,9 +948,9 @@ public:
* @param p Pressure (Pa) * @param p Pressure (Pa)
* @param tol Optional parameter setting the tolerance of the calculation. * @param tol Optional parameter setting the tolerance of the calculation.
* Important for some applications where numerical Jacobians * Important for some applications where numerical Jacobians
* are being calculated.Defaults to 1.0E-4 * are being calculated.
*/ */
virtual void setState_HP(doublereal h, doublereal p, doublereal tol = 1.e-4); virtual void setState_HP(double h, double p, double tol=1e-9);
//! Set the specific internal energy (J/kg) and specific volume (m^3/kg). //! Set the specific internal energy (J/kg) and specific volume (m^3/kg).
/*! /*!
@ -962,9 +962,9 @@ public:
* @param v specific volume (m^3/kg). * @param v specific volume (m^3/kg).
* @param tol Optional parameter setting the tolerance of the calculation. * @param tol Optional parameter setting the tolerance of the calculation.
* Important for some applications where numerical Jacobians * Important for some applications where numerical Jacobians
* are being calculated.Defaults to 1.0E-4 * are being calculated.
*/ */
virtual void setState_UV(doublereal u, doublereal v, doublereal tol = 1.e-4); virtual void setState_UV(double u, double v, double tol=1e-9);
//! Set the specific entropy (J/kg/K) and pressure (Pa). //! Set the specific entropy (J/kg/K) and pressure (Pa).
/*! /*!
@ -974,10 +974,10 @@ public:
* @param s specific entropy (J/kg/K) * @param s specific entropy (J/kg/K)
* @param p specific pressure (Pa). * @param p specific 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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_SP(doublereal s, doublereal p, doublereal tol = 1.e-4); virtual void setState_SP(double s, double p, double tol=1e-9);
//! Set the specific entropy (J/kg/K) and specific volume (m^3/kg). //! Set the specific entropy (J/kg/K) and specific volume (m^3/kg).
/*! /*!
@ -987,10 +987,10 @@ public:
* @param s specific entropy (J/kg/K) * @param s specific entropy (J/kg/K)
* @param v specific volume (m^3/kg). * @param v specific volume (m^3/kg).
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_SV(doublereal s, doublereal v, doublereal tol = 1.e-4); virtual void setState_SV(double s, double v, double tol=1e-9);
//! Set the specific entropy (J/kg/K) and temperature (K). //! Set the specific entropy (J/kg/K) and temperature (K).
/*! /*!
@ -1001,10 +1001,10 @@ public:
* @param s specific entropy (J/kg/K) * @param s specific entropy (J/kg/K)
* @param t temperature (K) * @param t temperature (K)
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_ST(double s, double t, double tol = 1.e-4) { virtual void setState_ST(double s, double t, double tol=1e-9) {
throw NotImplementedError("ThermoPhase::setState_ST"); throw NotImplementedError("ThermoPhase::setState_ST");
} }
@ -1017,10 +1017,10 @@ public:
* @param t temperature (K) * @param t temperature (K)
* @param v specific volume (m^3/kg) * @param v specific volume (m^3/kg)
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_TV(double t, double v, double tol = 1.e-4) { virtual void setState_TV(double t, double v, double tol=1e-9) {
throw NotImplementedError("ThermoPhase::setState_TV"); throw NotImplementedError("ThermoPhase::setState_TV");
} }
@ -1033,10 +1033,10 @@ public:
* @param p pressure (Pa) * @param p pressure (Pa)
* @param v specific volume (m^3/kg) * @param v specific volume (m^3/kg)
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_PV(double p, double v, double tol = 1.e-4) { virtual void setState_PV(double p, double v, double tol=1e-9) {
throw NotImplementedError("ThermoPhase::setState_PV"); throw NotImplementedError("ThermoPhase::setState_PV");
} }
@ -1049,10 +1049,10 @@ public:
* @param u specific internal energy (J/kg) * @param u specific internal energy (J/kg)
* @param p pressure (Pa) * @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_UP(double u, double p, double tol = 1.e-4) { virtual void setState_UP(double u, double p, double tol=1e-9) {
throw NotImplementedError("ThermoPhase::setState_UP"); throw NotImplementedError("ThermoPhase::setState_UP");
} }
@ -1065,10 +1065,10 @@ public:
* @param v specific volume (m^3/kg) * @param v specific volume (m^3/kg)
* @param h specific enthalpy (J/kg) * @param h specific enthalpy (J/kg)
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_VH(double v, double h, double tol = 1.e-4) { virtual void setState_VH(double v, double h, double tol=1e-9) {
throw NotImplementedError("ThermoPhase::setState_VH"); throw NotImplementedError("ThermoPhase::setState_VH");
} }
@ -1081,10 +1081,10 @@ public:
* @param t temperature (K) * @param t temperature (K)
* @param h specific enthalpy (J/kg) * @param h specific enthalpy (J/kg)
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_TH(double t, double h, double tol = 1.e-4) { virtual void setState_TH(double t, double h, double tol=1e-9) {
throw NotImplementedError("ThermoPhase::setState_TH"); throw NotImplementedError("ThermoPhase::setState_TH");
} }
@ -1097,10 +1097,10 @@ public:
* @param s specific entropy (J/kg/K) * @param s specific entropy (J/kg/K)
* @param h specific enthalpy (J/kg) * @param h specific enthalpy (J/kg)
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
*/ */
virtual void setState_SH(double s, double h, double tol = 1.e-4) { virtual void setState_SH(double s, double h, double tol=1e-9) {
throw NotImplementedError("ThermoPhase::setState_SH"); throw NotImplementedError("ThermoPhase::setState_SH");
} }
@ -1209,24 +1209,23 @@ private:
* @param h Specific enthalpy or internal energy (J/kg) * @param h Specific enthalpy or internal energy (J/kg)
* @param p Pressure (Pa) or specific volume (m^3/kg) * @param p Pressure (Pa) or specific volume (m^3/kg)
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
* @param doUV True if solving for UV, false for HP. * @param doUV True if solving for UV, false for HP.
*/ */
void setState_HPorUV(doublereal h, doublereal p, void setState_HPorUV(doublereal h, doublereal p,
doublereal tol = 1.e-4, bool doUV = false); doublereal tol=1e-9, bool doUV = false);
//! Carry out work in SP and SV calculations. //! Carry out work in SP and SV calculations.
/*! /*!
* @param s Specific entropy (J/kg) * @param s Specific entropy (J/kg)
* @param p Pressure (Pa) or specific volume (m^3/kg) * @param p Pressure (Pa) or specific volume (m^3/kg)
* @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. Important for some applications where * Important for some applications where numerical Jacobians
* numerical Jacobians are being calculated. * are being calculated.
* @param doSV True if solving for SV, false for SP. * @param doSV True if solving for SV, false for SP.
*/ */
void setState_SPorSV(doublereal s, doublereal p, void setState_SPorSV(double s, double p, double tol=1e-9, bool doSV = false);
doublereal tol = 1.e-4, bool doSV = false);
//! Helper function used by setState_HPorUV and setState_SPorSV. //! Helper function used by setState_HPorUV and setState_SPorSV.
//! Sets the temperature and (if set_p is true) the pressure. //! Sets the temperature and (if set_p is true) the pressure.

8
src/thermo/SingleSpeciesTP.cpp
View file

@ -215,7 +215,7 @@ void SingleSpeciesTP::setState_HP(doublereal h, doublereal p,
for (int n = 0; n < 50; n++) { for (int n = 0; n < 50; n++) {
dt = clip((h - enthalpy_mass())/cp_mass(), -100.0, 100.0); dt = clip((h - enthalpy_mass())/cp_mass(), -100.0, 100.0);
setState_TP(temperature() + dt, p); setState_TP(temperature() + dt, p);
if (fabs(dt) < tol) { if (fabs(dt / temperature()) < tol) {
return; return;
} }
} }
@ -234,7 +234,7 @@ void SingleSpeciesTP::setState_UV(doublereal u, doublereal v,
for (int n = 0; n < 50; n++) { for (int n = 0; n < 50; n++) {
dt = clip((u - intEnergy_mass())/cv_mass(), -100.0, 100.0); dt = clip((u - intEnergy_mass())/cv_mass(), -100.0, 100.0);
setTemperature(temperature() + dt); setTemperature(temperature() + dt);
if (fabs(dt) < tol) { if (fabs(dt / temperature()) < tol) {
return; return;
} }
} }
@ -250,7 +250,7 @@ void SingleSpeciesTP::setState_SP(doublereal s, doublereal p,
for (int n = 0; n < 50; n++) { for (int n = 0; n < 50; n++) {
dt = clip((s - entropy_mass())*temperature()/cp_mass(), -100.0, 100.0); dt = clip((s - entropy_mass())*temperature()/cp_mass(), -100.0, 100.0);
setState_TP(temperature() + dt, p); setState_TP(temperature() + dt, p);
if (fabs(dt) < tol) { if (fabs(dt / temperature()) < tol) {
return; return;
} }
} }
@ -269,7 +269,7 @@ void SingleSpeciesTP::setState_SV(doublereal s, doublereal v,
for (int n = 0; n < 50; n++) { for (int n = 0; n < 50; n++) {
dt = clip((s - entropy_mass())*temperature()/cv_mass(), -100.0, 100.0); dt = clip((s - entropy_mass())*temperature()/cv_mass(), -100.0, 100.0);
setTemperature(temperature() + dt); setTemperature(temperature() + dt);
if (fabs(dt) < tol) { if (fabs(dt / temperature()) < tol) {
return; return;
} }
} }

31
src/thermo/ThermoPhase.cpp
View file

@ -241,15 +241,14 @@ void ThermoPhase::setState_PY(doublereal p, doublereal* y)
setPressure(p); setPressure(p);
} }
void ThermoPhase::setState_HP(doublereal Htarget, doublereal p, void ThermoPhase::setState_HP(double Htarget, double p, double rtol)
doublereal dTtol)
{ {
setState_HPorUV(Htarget, p, dTtol, false); setState_HPorUV(Htarget, p, rtol, false);
} }
void ThermoPhase::setState_UV(doublereal u, doublereal v, doublereal dTtol) void ThermoPhase::setState_UV(double u, double v, double rtol)
{ {
setState_HPorUV(u, v, dTtol, true); setState_HPorUV(u, v, rtol, true);
} }
void ThermoPhase::setState_conditional_TP(doublereal t, doublereal p, bool set_p) void ThermoPhase::setState_conditional_TP(doublereal t, doublereal p, bool set_p)
@ -260,8 +259,8 @@ void ThermoPhase::setState_conditional_TP(doublereal t, doublereal p, bool set_p
} }
} }
void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p, void ThermoPhase::setState_HPorUV(double Htarget, double p,
doublereal dTtol, bool doUV) double rtol, bool doUV)
{ {
doublereal dt; doublereal dt;
doublereal v = 0.0; doublereal v = 0.0;
@ -411,7 +410,7 @@ void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
double acpd = std::max(fabs(cpd), 1.0E-5); double acpd = std::max(fabs(cpd), 1.0E-5);
double denom = std::max(fabs(Htarget), acpd * Tnew); double denom = std::max(fabs(Htarget), acpd * Tnew);
double HConvErr = fabs((Herr)/denom); double HConvErr = fabs((Herr)/denom);
if (HConvErr < 0.00001 * dTtol || fabs(dt/Tnew) < 0.00001 * dTtol) { if (HConvErr < rtol || fabs(dt/Tnew) < rtol) {
return; return;
} }
} }
@ -451,20 +450,18 @@ void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
} }
} }
void ThermoPhase::setState_SP(doublereal Starget, doublereal p, void ThermoPhase::setState_SP(double Starget, double p, double rtol)
doublereal dTtol)
{ {
setState_SPorSV(Starget, p, dTtol, false); setState_SPorSV(Starget, p, rtol, false);
} }
void ThermoPhase::setState_SV(doublereal Starget, doublereal v, void ThermoPhase::setState_SV(double Starget, double v, double rtol)
doublereal dTtol)
{ {
setState_SPorSV(Starget, v, dTtol, true); setState_SPorSV(Starget, v, rtol, true);
} }
void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p, void ThermoPhase::setState_SPorSV(double Starget, double p,
doublereal dTtol, bool doSV) double rtol, bool doSV)
{ {
doublereal v = 0.0; doublereal v = 0.0;
doublereal dt; doublereal dt;
@ -597,7 +594,7 @@ void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p,
double acpd = std::max(fabs(cpd), 1.0E-5); double acpd = std::max(fabs(cpd), 1.0E-5);
double denom = std::max(fabs(Starget), acpd * Tnew); double denom = std::max(fabs(Starget), acpd * Tnew);
double SConvErr = fabs((Serr * Tnew)/denom); double SConvErr = fabs((Serr * Tnew)/denom);
if (SConvErr < 0.00001 * dTtol || fabs(dt/Tnew) < 0.00001 * dTtol) { if (SConvErr < rtol || fabs(dt/Tnew) < rtol) {
return; return;
} }
} }

2
src/zeroD/ConstPressureReactor.cpp
View file

@ -57,7 +57,7 @@ void ConstPressureReactor::updateState(doublereal* y)
m_mass = y[0]; m_mass = y[0];
m_thermo->setMassFractions_NoNorm(y+2); m_thermo->setMassFractions_NoNorm(y+2);
if (m_energy) { if (m_energy) {
m_thermo->setState_HP(y[1]/m_mass, m_pressure, 1.0e-4); m_thermo->setState_HP(y[1]/m_mass, m_pressure);
} else { } else {
m_thermo->setPressure(m_pressure); m_thermo->setPressure(m_pressure);
} }