From 0dcdf7ef5dda49625854260869502c3b880529fe Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 20 Aug 2014 21:44:13 +0200 Subject: [PATCH] Implemented PQ for pseudo-pure to close https://github.com/CoolProp/CoolProp/issues/92 * As a consequence, made update function take CoolProp::input_pairs enum as first input (thanks @jowr for the idea) * Decoupled the pre_update and post_update functions in order to allow for more creative update functions using guess values, etc. Signed-off-by: Ian Bell --- include/AbstractState.h | 2 +- include/DataStructures.h | 4 +- src/Backends/Helmholtz/FlashRoutines.cpp | 18 +++++++- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 45 +++++++++++++++---- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 9 +++- .../Incompressible/IncompressibleBackend.cpp | 2 +- .../Incompressible/IncompressibleBackend.h | 2 +- .../REFPROP/REFPROPMixtureBackend.cpp | 2 +- src/Backends/REFPROP/REFPROPMixtureBackend.h | 2 +- src/Backends/Tabular/TabularBackends.h | 2 +- src/CoolProp.cpp | 2 +- src/SpeedTest.cpp | 2 +- src/Tests/CoolProp-Tests.cpp | 15 ++++--- 13 files changed, 78 insertions(+), 29 deletions(-) diff --git a/include/AbstractState.h b/include/AbstractState.h index e74814a4..db7ac522 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -274,7 +274,7 @@ public: virtual bool using_mass_fractions(void) = 0; virtual bool using_volu_fractions(void) = 0; - virtual void update(long input_pair, double Value1, double Value2) = 0; + virtual void update(CoolProp::input_pairs input_pair, double Value1, double Value2) = 0; virtual void set_mole_fractions(const std::vector &mole_fractions) = 0; virtual void set_mass_fractions(const std::vector &mass_fractions) = 0; virtual void set_volu_fractions(const std::vector &mass_fractions){throw NotImplementedError("Volume composition has not been implemented.");} diff --git a/include/DataStructures.h b/include/DataStructures.h index ac2751bd..5325c645 100644 --- a/include/DataStructures.h +++ b/include/DataStructures.h @@ -152,9 +152,9 @@ inline bool match_pair(long key1, long key2, long x1, long x2, bool &swap) swap = !(key1 == x1); return ((key1 == x1 && key2 == x2) || (key2 == x1 && key1 == x2)); }; -template long generate_update_pair(long key1, T value1, long key2, T value2, T &out1, T&out2) +template CoolProp::input_pairs generate_update_pair(long key1, T value1, long key2, T value2, T &out1, T&out2) { - long pair; + CoolProp::input_pairs pair; bool swap; if (match_pair(key1, key2, iQ, iT, swap)){ diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index e487b610..89db0818 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -140,8 +140,22 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) if (HEOS.is_pure_or_pseudopure) { if (HEOS.components[0]->pEOS->pseudo_pure){ - // It is a psedo-pure mixture - throw NotImplementedError("PQ_flash not implemented for pseudo-pure fluids yet"); + // It is a pseudo-pure mixture + + HEOS._TLanc = HEOS.components[0]->ancillaries.pL.invert(HEOS._p); + HEOS._TVanc = HEOS.components[0]->ancillaries.pV.invert(HEOS._p); + // Get guesses for the ancillaries for density + long double rhoL = HEOS.components[0]->ancillaries.rhoL.evaluate(HEOS._TLanc); + long double rhoV = HEOS.components[0]->ancillaries.rhoV.evaluate(HEOS._TVanc); + // Solve for the density + HEOS.SatL->update_TP_guessrho(HEOS._TLanc, HEOS._p, rhoL); + HEOS.SatV->update_TP_guessrho(HEOS._TVanc, HEOS._p, rhoV); + + // Load the outputs + HEOS._phase = iphase_twophase; + HEOS._p = HEOS._Q*HEOS.SatV->p() + (1 - HEOS._Q)*HEOS.SatL->p(); + HEOS._T = HEOS._Q*HEOS.SatV->T() + (1 - HEOS._Q)*HEOS.SatL->T(); + HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar()); } else{ // Critical point for pure fluids, slightly different for pseudo-pure, very different for mixtures diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index c5f1c82b..53f11bf4 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -535,13 +535,23 @@ long double HelmholtzEOSMixtureBackend::calc_pmax(void) return summer; } -void HelmholtzEOSMixtureBackend::update_TP_guessrho(long double T, long double p, long double rho_guess) +void HelmholtzEOSMixtureBackend::update_TP_guessrho(long double T, long double p, long double rhomolar_guess) { - double rho = solver_rho_Tp(T, p, rho_guess); - update(DmolarT_INPUTS, rho, T); + CoolProp::input_pairs pair = PT_INPUTS; + // Set up the state + pre_update(pair, p, T); + + // Do the flash call + double rhomolar = solver_rho_Tp(T, p, rhomolar_guess); + + // Update the class with the new calculated density + update(DmolarT_INPUTS, rhomolar, T); + + // Cleanup + post_update(); } -void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double &value1, double &value2) +void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(CoolProp::input_pairs &input_pair, long double &value1, long double &value2) { // Check if a mass based input, convert it to molar units @@ -588,24 +598,33 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double & return; } } -void HelmholtzEOSMixtureBackend::update(long input_pair, double value1, double value2 ) + +void HelmholtzEOSMixtureBackend::pre_update(CoolProp::input_pairs &input_pair, long double &value1, long double &value2 ) { - if (get_debug_level() > 10){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 the state clear(); if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) { throw ValueError("Mole fractions must be set"); } - + + // If the inputs are in mass units, convert them to molar units mass_to_molar_inputs(input_pair, value1, value2); // Set the mole-fraction weighted gas constant for the mixture // (or the pure/pseudo-pure fluid) if it hasn't been set yet gas_constant(); - // Reducing state + // Calculate and cache the reducing state calc_reducing_state(); +} + +void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double value1, double value2 ) +{ + if (get_debug_level() > 10){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;} + + long double ld_value1 = value1, ld_value2 = value2; + pre_update(input_pair, ld_value1, ld_value2); switch(input_pair) { @@ -640,6 +659,13 @@ void HelmholtzEOSMixtureBackend::update(long input_pair, double value1, double v default: throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str())); } + + post_update(); + +} + +void HelmholtzEOSMixtureBackend::post_update() +{ // Check the values that must always be set //if (_p < 0){ throw ValueError("p is less than zero");} if (!ValidNumber(_p)){ throw ValueError("p is not a valid number");} @@ -1960,6 +1986,7 @@ SimpleState HelmholtzEOSMixtureBackend::calc_reducing_state_nocache(const std::v } void HelmholtzEOSMixtureBackend::calc_reducing_state(void) { + /// \todo set critical independently _reducing = calc_reducing_state_nocache(mole_fractions); _crit = _reducing; } diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index 86c1947f..21dfb6d4 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -15,6 +15,10 @@ namespace CoolProp { class FlashRoutines; class HelmholtzEOSMixtureBackend : public AbstractState { + +private: + void pre_update(CoolProp::input_pairs &input_pair, long double &value1, long double &value2 ); + void post_update(); protected: std::vector components; ///< The components that are in use @@ -26,6 +30,7 @@ protected: SimpleState _crit; phases imposed_phase_index; int N; ///< Number of components + public: HelmholtzEOSMixtureBackend(){imposed_phase_index = iphase_not_imposed; _phase = iphase_unknown;}; HelmholtzEOSMixtureBackend(std::vector components, bool generate_SatL_and_SatV = true); @@ -58,7 +63,7 @@ public: void resize(unsigned int N); shared_ptr SatL, SatV; ///< - void update(long input_pair, double value1, double value2); + void update(CoolProp::input_pairs input_pair, double value1, double value2); void update_TP_guessrho(long double T, long double p, long double rho_guess); @@ -209,7 +214,7 @@ public: void update_states(); - void mass_to_molar_inputs(long &input_pair, double &value1, double &value2); + void mass_to_molar_inputs(CoolProp::input_pairs &input_pair, long double &value1, long double &value2); // *************************************************************** // *************************************************************** diff --git a/src/Backends/Incompressible/IncompressibleBackend.cpp b/src/Backends/Incompressible/IncompressibleBackend.cpp index 8b3a4bfa..36cc658d 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.cpp +++ b/src/Backends/Incompressible/IncompressibleBackend.cpp @@ -37,7 +37,7 @@ IncompressibleBackend::IncompressibleBackend(const std::vector &com throw NotImplementedError("Mixture-style constructor is not implemented yet for incompressible fluids"); } -void IncompressibleBackend::update(long input_pair, double value1, double value2) { +void IncompressibleBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) { //if (mass_fractions.empty()){ // throw ValueError("mass fractions have not been set"); //} diff --git a/src/Backends/Incompressible/IncompressibleBackend.h b/src/Backends/Incompressible/IncompressibleBackend.h index 43af6d64..7f23f54f 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.h +++ b/src/Backends/Incompressible/IncompressibleBackend.h @@ -46,7 +46,7 @@ public: @param value1 First input value @param value2 Second input value */ - void update(long input_pair, double value1, double value2); + void update(CoolProp::input_pairs input_pair, double value1, double value2); /// Set the mole fractions /** diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp index ccd85a32..e791c0d0 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp @@ -736,7 +736,7 @@ void REFPROPMixtureBackend::calc_phase_envelope(const std::string &type) if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } } -void REFPROPMixtureBackend::update(long input_pair, double value1, double value2) +void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) { double rho_mol_L=_HUGE, rhoLmol_L=_HUGE, rhoVmol_L=_HUGE, hmol=_HUGE,emol=_HUGE,smol=_HUGE,cvmol=_HUGE,cpmol=_HUGE, diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.h b/src/Backends/REFPROP/REFPROPMixtureBackend.h index 35d5a800..3014d54b 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.h +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.h @@ -46,7 +46,7 @@ public: @param value1 First input value @param value2 Second input value */ - void update(long input_pair, + void update(CoolProp::input_pairs, double value1, double value2 ); diff --git a/src/Backends/Tabular/TabularBackends.h b/src/Backends/Tabular/TabularBackends.h index db437526..3b54363b 100644 --- a/src/Backends/Tabular/TabularBackends.h +++ b/src/Backends/Tabular/TabularBackends.h @@ -43,7 +43,7 @@ class GriddedTableBackend : public AbstractState bool using_mass_fractions(void){return false;} bool using_volu_fractions(void){return false;} - void update(long input_pair, double Value1, double Value2){}; + void update(CoolProp::input_pairs input_pair, double Value1, double Value2){}; void set_mole_fractions(const std::vector &mole_fractions){}; void set_mass_fractions(const std::vector &mass_fractions){}; diff --git a/src/CoolProp.cpp b/src/CoolProp.cpp index ca4c6851..c77c6e88 100644 --- a/src/CoolProp.cpp +++ b/src/CoolProp.cpp @@ -386,7 +386,7 @@ double _PropsSI(const std::string &Output, const std::string &Name1, double Prop } // Obtain the input pair - long pair = generate_update_pair(iName1, Prop1, iName2, Prop2, x1, x2); + CoolProp::input_pairs pair = generate_update_pair(iName1, Prop1, iName2, Prop2, x1, x2); // Update the state State->update(pair, x1, x2); diff --git a/src/SpeedTest.cpp b/src/SpeedTest.cpp index 2d8208c8..1f0e145a 100644 --- a/src/SpeedTest.cpp +++ b/src/SpeedTest.cpp @@ -8,7 +8,7 @@ namespace CoolProp{ -void compare_REFPROP_and_CoolProp(std::string fluid, int inputs, double val1, double val2, std::size_t N, double d1, double d2) +void compare_REFPROP_and_CoolProp(std::string fluid, CoolProp::input_pairs inputs, double val1, double val2, std::size_t N, double d1, double d2) { time_t t1,t2; diff --git a/src/Tests/CoolProp-Tests.cpp b/src/Tests/CoolProp-Tests.cpp index 238f3cc5..8302ecbc 100644 --- a/src/Tests/CoolProp-Tests.cpp +++ b/src/Tests/CoolProp-Tests.cpp @@ -221,7 +221,7 @@ class TransportValidationFixture protected: long double actual, x1, x2; shared_ptr pState; - int pair; + CoolProp::input_pairs pair; public: TransportValidationFixture(){ } ~TransportValidationFixture(){ } @@ -232,7 +232,7 @@ public: double o1, o2; long iin1 = CoolProp::get_parameter_index(in1); long iin2 = CoolProp::get_parameter_index(in2); - long pair = CoolProp::generate_update_pair(iin1, v1, iin2, v2, o1, o2); + CoolProp::input_pairs pair = CoolProp::generate_update_pair(iin1, v1, iin2, v2, o1, o2); pState->update(pair, o1, o2); } void get_value(long key) @@ -484,7 +484,7 @@ TEST_CASE_METHOD(TransportValidationFixture, "Compare thermal conductivities aga }; /* namespace TransportValidation */ -static int inputs[] = { +static CoolProp::input_pairs inputs[] = { CoolProp::DmolarT_INPUTS, //CoolProp::SmolarT_INPUTS, //CoolProp::HmolarT_INPUTS, @@ -511,14 +511,14 @@ class ConsistencyFixture protected: long double hmolar, pmolar, smolar, umolar, rhomolar, T, p, x1, x2; shared_ptr pState; - int pair; + CoolProp::input_pairs pair; public: ConsistencyFixture(){} ~ConsistencyFixture(){} void set_backend(std::string backend, std::string fluid_name){ pState.reset(CoolProp::AbstractState::factory(backend, fluid_name)); } - void set_pair(int pair){ + void set_pair(CoolProp::input_pairs pair){ this->pair = pair; } void set_TP(long double T, long double p) @@ -569,6 +569,9 @@ public: x1 = hmolar; x2 = smolar; break; case CoolProp::SmolarUmolar_INPUTS: x1 = smolar; x2 = umolar; break; + + default: + throw CoolProp::ValueError(); } } void single_phase_consistency_check() @@ -595,7 +598,7 @@ TEST_CASE_METHOD(ConsistencyFixture, "Test all input pairs for CO2 using all val for (int i = 0; i < inputsN; ++i) { - int pair = inputs[i]; + CoolProp::input_pairs pair = inputs[i]; std::string pair_desc = CoolProp::get_input_pair_short_desc(pair); set_pair(pair); CAPTURE(pair_desc);