diff --git a/include/AbstractState.h b/include/AbstractState.h index 7ddd4f7b..b266b459 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -39,7 +39,7 @@ protected: /// Some administrative variables long _fluid_type; - long _phase; ///< The key for the phase from CoolProp::phases enum + phases _phase; ///< The key for the phase from CoolProp::phases enum bool _forceSinglePhase, _forceTwoPhase; bool isSupercriticalPhase(void){ @@ -185,7 +185,9 @@ protected: virtual long double calc_physical_hazard(void){throw NotImplementedError("calc_physical_hazard is not implemented for this backend");}; /// Calculate the first partial derivative for the desired derivative - virtual long double calc_first_partial_deriv(int Of, int Wrt, int Constant){throw NotImplementedError("calc_first_partial_deriv is not implemented for this backend");}; + virtual long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant){throw NotImplementedError("calc_first_partial_deriv is not implemented for this backend");}; + /// Calculate the second partial derivative using the given backend + virtual long double calc_second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Of2, parameters Constant2){throw NotImplementedError("calc_second_partial_deriv is not implemented for this backend");}; /// Using this backend, calculate the reduced density (rho/rhoc) virtual long double calc_reduced_density(void){throw NotImplementedError("calc_reduced_density is not implemented for this backend");}; @@ -232,7 +234,7 @@ protected: virtual long double calc_melting_line(int param, int given, long double value){throw NotImplementedError("This backend does not implement calc_melting_line function");}; - virtual int calc_phase(void){throw NotImplementedError("This backend does not implement calc_phase function");}; + virtual phases calc_phase(void){throw NotImplementedError("This backend does not implement calc_phase function");}; /// Using this backend, calculate a phase given by the state string /// @param state A string that describes the state desired, one of "hs_anchor", "critical"/"crit", "reducing" @@ -291,8 +293,31 @@ public: double keyed_output(int key); double trivial_keyed_output(int key); - - long double first_partial_deriv(int Of, int Wrt, int Constant){return calc_first_partial_deriv(Of,Wrt,Constant);}; + + /** \brief The first partial derivative in homogeneous phases + * + * \f[ \left(\frac{\partial A}{\partial B}\right)_C = \frac{\left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta}{\left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta} = \frac{N}{D}\f] + */ + long double first_partial_deriv(parameters Of, parameters Wrt, parameters Constant){return calc_first_partial_deriv(Of, Wrt, Constant);}; + + /** \brief The second partial derivative in homogeneous phases + * + * \sa \ref CoolProp::AbstractState::first_partial_deriv + * + * \f[ + * \frac{\partial}{\partial D}\left(\left(\frac{\partial A}{\partial B}\right)_C\right)_E = \frac{\frac{\partial}{\partial \tau}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\delta\left(\frac{\partial E}{\partial \delta}\right)_\tau-\frac{\partial}{\partial \delta}\left(\left(\frac{\partial A}{\partial B}\right)_C\right)_\tau\left(\frac{\partial E}{\partial \tau}\right)_\delta}{\left(\frac{\partial D}{\partial \tau}\right)_\delta\left(\frac{\partial E}{\partial \delta}\right)_\tau-\left(\frac{\partial D}{\partial \delta}\right)_\tau\left(\frac{\partial E}{\partial \tau}\right)_\delta} + * \f] + * + * which can be expressed in parts as + * + * \f[\left(\frac{\partial N}{\partial \delta}\right)_{\tau} = \left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta^2}\right)_{\tau}+\left(\frac{\partial^2 A}{\partial \tau\partial\delta}\right)\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau\partial\delta}\right)-\left(\frac{\partial^2 A}{\partial \delta^2}\right)_{\tau}\left(\frac{\partial C}{\partial \tau}\right)_\delta\f] + * \f[\left(\frac{\partial D}{\partial \delta}\right)_{\tau} = \left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta^2}\right)_{\tau}+\left(\frac{\partial^2 B}{\partial \tau\partial\delta}\right)\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau\partial\delta}\right)-\left(\frac{\partial^2 B}{\partial \delta^2}\right)_{\tau}\left(\frac{\partial C}{\partial \tau}\right)_\delta\f] + * \f[\left(\frac{\partial N}{\partial \tau}\right)_{\delta} = \left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta\partial\tau}\right)+\left(\frac{\partial^2 A}{\partial \tau^2}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau^2}\right)_\delta-\left(\frac{\partial^2 A}{\partial \delta\partial\tau}\right)\left(\frac{\partial C}{\partial \tau}\right)_\delta\f] + * \f[\left(\frac{\partial D}{\partial \tau}\right)_{\delta} = \left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta\partial\tau}\right)+\left(\frac{\partial^2 B}{\partial \tau^2}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau^2}\right)_\delta-\left(\frac{\partial^2 B}{\partial \delta\partial\tau}\right)\left(\frac{\partial C}{\partial \tau}\right)_\delta\f] + * \f[\frac{\partial}{\partial \tau}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\delta = \frac{D\left(\frac{\partial N}{\partial \tau}\right)_{\delta}-N\left(\frac{\partial D}{\partial \tau}\right)_{\delta}}{D^2}\f] + * \f[\frac{\partial}{\partial \delta}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\tau = \frac{D\left(\frac{\partial N}{\partial \delta}\right)_{\tau}-N\left(\frac{\partial D}{\partial \delta}\right)_{\tau}}{D^2}\f] + */ + long double second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Of2, parameters Constant2){return calc_second_partial_deriv(Of1,Wrt1,Constant1,Of2,Constant2);}; // Limits double Tmin(void); @@ -300,7 +325,7 @@ public: double pmax(void); double Ttriple(void); - int phase(void){return calc_phase();}; + phases phase(void){return calc_phase();}; /// Return the critical temperature in K /** diff --git a/include/DataStructures.h b/include/DataStructures.h index 9f79be52..81129625 100644 --- a/include/DataStructures.h +++ b/include/DataStructures.h @@ -92,7 +92,8 @@ enum phases{iphase_liquid, ///< Subcritical liquid iphase_supercritical_liquid, ///< Supercritical liquid (p > pc, T < Tc) iphase_gas, ///< Subcritical gas iphase_twophase, ///< Twophase - iphase_unknown}; + iphase_unknown, ///< Unknown phase + iphase_not_imposed}; ///< Phase is not imposed /// These are unit types for the fluid enum fluid_types{FLUID_TYPE_PURE, FLUID_TYPE_PSEUDOPURE, FLUID_TYPE_REFPROP, FLUID_TYPE_INCOMPRESSIBLE_LIQUID, FLUID_TYPE_INCOMPRESSIBLE_SOLUTION, FLUID_TYPE_UNDEFINED}; diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 7d0e15c2..241c8b37 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -12,6 +12,7 @@ #include "Backends/Helmholtz/HelmholtzEOSBackend.h" #include "Backends/Incompressible/IncompressibleBackend.h" #include "Backends/Helmholtz/Fluids/FluidLibrary.h" +#include "Backends/Tabular/TabularBackends.h" namespace CoolProp { @@ -46,9 +47,11 @@ AbstractState * AbstractState::factory(const std::string &backend, const std::st { return new IncompressibleBackend(fluid_string); } - else if (!backend.compare("BRINE")) + else if (backend.find("TTSE&") == 0) { - throw ValueError("BRINE backend not yet implemented"); + // Will throw if there is a problem with this backend + shared_ptr AS(factory(backend.substr(5), fluid_string)); + return new TTSEBackend(*AS.get()); } else if (!backend.compare("TREND")) { diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index 250033db..0835dc1f 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -6,7 +6,7 @@ namespace CoolProp{ void FlashRoutines::PT_flash(HelmholtzEOSMixtureBackend &HEOS) { - if (HEOS.imposed_phase_index == -1) // If no phase index is imposed (see set_components function) + if (HEOS.imposed_phase_index == iphase_not_imposed) // If no phase index is imposed (see set_components function) { // Find the phase, while updating all internal variables possible HEOS.T_phase_determination_pure_or_pseudopure(iP, HEOS._p); @@ -25,6 +25,7 @@ void FlashRoutines::PT_flash(HelmholtzEOSMixtureBackend &HEOS) void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS) { + long double T = HEOS._T; if (HEOS.is_pure_or_pseudopure) { // The maximum possible saturation temperature @@ -67,8 +68,13 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS) } } catch(std::exception &){ - // We may need to polish the solution at low pressure - SaturationSolvers::saturation_T_pure_1D_P(&HEOS, HEOS._T, options); + try{ + // We may need to polish the solution at low pressure + SaturationSolvers::saturation_T_pure_1D_P(&HEOS, T, options); + } + catch(std::exception &){ + SaturationSolvers::saturation_critical(&HEOS, iT, T); + } } // Load the outputs HEOS._phase = iphase_twophase; @@ -263,7 +269,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other) std::string errstring; - if (HEOS.imposed_phase_index > -1) + if (HEOS.imposed_phase_index != iphase_not_imposed) { // Use the phase defined by the imposed phase HEOS._phase = HEOS.imposed_phase_index; @@ -467,7 +473,7 @@ void FlashRoutines::HSU_P_flash_singlephase_Newton(HelmholtzEOSMixtureBackend &H break; } default: - break; + throw ValueError("other variable in HSU_P_flash_singlephase_Newton is invalid"); } //First index is the row, second index is the column @@ -563,7 +569,7 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other) { bool saturation_called = false; long double value; - if (HEOS.imposed_phase_index > -1) + if (HEOS.imposed_phase_index != iphase_not_imposed) { // Use the phase defined by the imposed phase HEOS._phase = HEOS.imposed_phase_index; @@ -640,7 +646,7 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other) } void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, int other) { - if (HEOS.imposed_phase_index > -1) + if (HEOS.imposed_phase_index != iphase_not_imposed) { // Use the phase defined by the imposed phase HEOS._phase = HEOS.imposed_phase_index; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 0282bf8f..e6a360cb 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -74,7 +74,7 @@ void HelmholtzEOSMixtureBackend::set_components(std::vector comp set_excess_term(); } - imposed_phase_index = -1; + imposed_phase_index = iphase_not_imposed; // Top-level class can hold copies of the base saturation classes, // saturation classes cannot hold copies of the saturation classes @@ -581,7 +581,6 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double & case DmassSmass_INPUTS: input_pair = DmolarSmolar_INPUTS; value1 /= mm; value2 *= mm; break; case DmassUmass_INPUTS: input_pair = DmolarUmolar_INPUTS; value1 /= mm; value2 *= mm; break; } - } default: return; @@ -589,6 +588,8 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double & } void HelmholtzEOSMixtureBackend::update(long input_pair, double value1, double value2 ) { + if (get_debug_level() > 0){std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)",__FILE__,__LINE__, input_pair, get_input_pair_short_desc(input_pair).c_str(), value1, value2) << std::endl;} + clear(); if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) { @@ -1305,7 +1306,7 @@ long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv_nocache(long do return (dOf_dtau*dConstant_ddelta-dOf_ddelta*dConstant_dtau)/(dWrt_dtau*dConstant_ddelta-dWrt_ddelta*dConstant_dtau); } -long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv(int Of, int Wrt, int Constant) +long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) { return calc_first_partial_deriv_nocache(_T, _rhomolar, Of, Wrt, Constant); } @@ -1465,7 +1466,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do } long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double p, long double rhomolar_guess) { - int phase; + phases phase; // Define the residual to be driven to zero class solver_TP_resid : public FuncWrapper1D @@ -1495,10 +1496,14 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double solver_TP_resid resid(this,T,p); std::string errstring; - if (imposed_phase_index > -1) + // Check if the phase is imposed + if (imposed_phase_index != iphase_not_imposed) + // Use the imposed phase index phase = imposed_phase_index; else + // Use the phase index in the class phase = _phase; + if (rhomolar_guess < 0) // Not provided { // Calculate a guess value using SRK equation of state diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index 78de62d8..3c8d0792 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -24,10 +24,10 @@ protected: lnK; ///< The natural logarithms of the K factors of the components SimpleState _crit; - int imposed_phase_index; + phases imposed_phase_index; int N; ///< Number of components public: - HelmholtzEOSMixtureBackend(){imposed_phase_index = -1; _phase = iphase_unknown;}; + HelmholtzEOSMixtureBackend(){imposed_phase_index = iphase_not_imposed; _phase = iphase_unknown;}; HelmholtzEOSMixtureBackend(std::vector components, bool generate_SatL_and_SatV = true); HelmholtzEOSMixtureBackend(std::vector &component_names, bool generate_SatL_and_SatV = true); virtual ~HelmholtzEOSMixtureBackend(){}; @@ -47,7 +47,7 @@ public: bool has_melting_line(){ return is_pure_or_pseudopure && components[0]->ancillaries.melting_line.enabled();}; long double calc_melting_line(int param, int given, long double value); - int calc_phase(void){return _phase;}; + phases calc_phase(void){return _phase;}; const CoolProp::SimpleState &calc_state(const std::string &state); @@ -73,7 +73,7 @@ public: \brief Specify the phase - this phase will always be used in calculations @param phase_index The index from CoolProp::phases */ - void specify_phase(int phase_index){imposed_phase_index = phase_index; _phase = phase_index;}; + void specify_phase(phases phase_index){imposed_phase_index = phase_index; _phase = phase_index;}; void set_reducing_function(); void set_excess_term(); @@ -193,7 +193,7 @@ public: \left(\frac{\partial A}{\partial B}\right)_C = \frac{\left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta}{\left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta} \f] */ - long double calc_first_partial_deriv(int Of, int Wrt, int Constant); + long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant); /** This version doesn't use any cached values diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 2764450e..f3c7cc84 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -4,8 +4,77 @@ #include "MatrixMath.h" namespace CoolProp { + +void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend *HEOS, parameters ykey, long double y){ + + class inner_resid : public FuncWrapper1D{ + public: + HelmholtzEOSMixtureBackend *HEOS; + long double T, desired_p, rhomolar_liq, calc_p, rhomolar_crit; + + inner_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double desired_p) + : HEOS(HEOS), T(T), desired_p(desired_p){}; + double call(double rhomolar_liq){ + this->rhomolar_liq = rhomolar_liq; + HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T); + calc_p = HEOS->SatL->p(); + std::cout << format("inner p: %0.16Lg; res: %0.16Lg", calc_p, calc_p - desired_p) << std::endl; + return calc_p - desired_p; + } + }; + + // Define the outer residual to be driven to zero - this is the equality of + // Gibbs function for both co-existing phases + class outer_resid : public FuncWrapper1D + { + public: - void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options){ + HelmholtzEOSMixtureBackend *HEOS; + parameters ykey; + long double y; + long double r, T, rhomolar_liq, rhomolar_vap, value, p, gL, gV, rhomolar_crit; + int other; + + outer_resid(HelmholtzEOSMixtureBackend *HEOS, CoolProp::parameters ykey, long double y) + : HEOS(HEOS), ykey(ykey), y(y){ + rhomolar_crit = HEOS->rhomolar_critical(); + }; + double call(double rhomolar_vap){ + this->y = y; + + // Calculate the other variable (T->p or p->T) for given vapor density + if (ykey == iT){ + T = y; + HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, y); + this->p = HEOS->SatV->p(); + std::cout << format("outer p: %0.16Lg",this->p) << std::endl; + inner_resid inner(HEOS, T, p); + std::string errstr2; + rhomolar_liq = Brent(inner, rhomolar_crit*1.5, rhomolar_crit*(1+1e-8), LDBL_EPSILON, 1e-10, 100, errstr2); + } + HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T); + HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, T); + + // Calculate the Gibbs functions for liquid and vapor + gL = HEOS->SatL->gibbsmolar(); + gV = HEOS->SatV->gibbsmolar(); + + // Residual is difference in Gibbs function + r = gL - gV; + + return this->p; + }; + }; + outer_resid resid(HEOS, iT, y); + + double rhomolar_crit = HEOS->rhomolar_critical(); + + std::string errstr; + Brent(&resid, rhomolar_crit*(1-1e-8), rhomolar_crit*0.5, DBL_EPSILON, 1e-9, 20, errstr); +} + +void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options) +{ // Define the residual to be driven to zero class solver_resid : public FuncWrapper1D diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index f690fa88..05b4d06c 100644 --- a/src/Backends/Helmholtz/VLERoutines.h +++ b/src/Backends/Helmholtz/VLERoutines.h @@ -15,8 +15,8 @@ namespace SaturationSolvers }; struct saturation_T_pure_options{ bool use_guesses; ///< true to start off at the values specified by rhoL, rhoV - long double omega, rhoL, rhoV, pL, pV, p; - saturation_T_pure_options(){omega = _HUGE; rhoV = _HUGE; rhoL = _HUGE; rhoL = _HUGE; pV = _HUGE, pL = _HUGE;} + long double omega, rhoL, rhoV, pL, pV, p, T; + saturation_T_pure_options(){omega = _HUGE; rhoV = _HUGE; rhoL = _HUGE; rhoL = _HUGE; pV = _HUGE, pL = _HUGE; T = _HUGE;} }; struct saturation_D_pure_options{ @@ -89,6 +89,21 @@ namespace SaturationSolvers */ void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options); + /* \brief A robust but slow solver in the very-near-critical region + * + * This solver operates in the following fashion: + * 1. Using a bounded interval for rho'':[rhoc, rhoc-??], guess a value for rho'' + * 2. For guessed value of rho'' and given value of T, calculate p + * 3. Using a Brent solver on the other co-existing phase (rho'), calculate the (bounded) value of rho' that yields the same pressure + * 4. Use another outer Brent solver on rho'' to enforce the same Gibbs function between liquid and vapor + * 5. Fin. + * + * @param HEOS The Helmholtz EOS backend instance to be used + * @param ykey The CoolProp::parameters key to be imposed - one of iT or iP + * @param y The value for the imposed variable + */ + void saturation_critical(HelmholtzEOSMixtureBackend *HEOS, CoolProp::parameters ykey, long double y); + void successive_substitution(HelmholtzEOSMixtureBackend &HEOS, const long double beta, long double T, diff --git a/src/Tests/CoolProp-Tests.cpp b/src/Tests/CoolProp-Tests.cpp index 97262225..a56558b3 100644 --- a/src/Tests/CoolProp-Tests.cpp +++ b/src/Tests/CoolProp-Tests.cpp @@ -793,4 +793,45 @@ TEST_CASE("Tests for solvers in P,Y flash using Water", "[flash],[PH],[PS],[PU]" } } +//TEST_CASE("Test that states agree with CoolProp", "[states]") +//{ +// std::vector fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),','); +// for (std::size_t i = 0; i < fluids.size(); ++i) +// { +// shared_ptr AS(CoolProp::AbstractState::factory("HEOS",fluids[i])); +// +// CoolProp::SimpleState &triple_liquid = AS->get_state("triple_liquid"); +// CoolProp::SimpleState &triple_vapor = AS->get_state("triple_vapor"); +// +// // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU +// std::ostringstream ss1; +// ss1 << "Check state for triple_liquid for " << fluids[i]; +// SECTION(ss1.str(),"") +// { +// std::string note = "The enthalpy and entropy are hardcoded in the fluid JSON files. They MUST agree with the values calculated by the EOS"; +// AS->update(CoolProp::DmolarT_INPUTS, hs_anchor.rhomolar, hs_anchor.T); +// CAPTURE(hs_anchor.hmolar); +// CAPTURE(hs_anchor.smolar); +// double EOS_hmolar = AS->hmolar(); +// double EOS_smolar = AS->smolar(); +// CHECK( std::abs(EOS_hmolar - hs_anchor.hmolar) < 1e-3); +// CHECK( std::abs(EOS_smolar - hs_anchor.smolar) < 1e-3); +// } +// +// // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU +// std::ostringstream ss2; +// ss2 << "Check state for triple_vapor for " << fluids[i]; +// SECTION(ss2.str(),"") +// { +// std::string note = "The enthalpy and entropy are hardcoded in the fluid JSON files. They MUST agree with the values calculated by the EOS"; +// AS->update(CoolProp::DmolarT_INPUTS, hs_anchor.rhomolar, hs_anchor.T); +// CAPTURE(hs_anchor.hmolar); +// CAPTURE(hs_anchor.smolar); +// double EOS_hmolar = AS->hmolar(); +// double EOS_smolar = AS->smolar(); +// CHECK( std::abs(EOS_hmolar - hs_anchor.hmolar) < 1e-3); +// CHECK( std::abs(EOS_smolar - hs_anchor.smolar) < 1e-3); +// } +// } +//} #endif