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 <ian.h.bell@gmail.com>

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<long double> &mole_fractions) = 0;
     virtual void set_mass_fractions(const std::vector<long double> &mass_fractions) = 0;
     virtual void set_volu_fractions(const std::vector<long double> &mass_fractions){throw NotImplementedError("Volume composition has not been implemented.");}

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<class T> long generate_update_pair(long key1, T value1, long key2, T value2, T &out1, T&out2)
+template<class T> 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)){

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

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;
 }

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<CoolPropFluid*> 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<CoolPropFluid*> components, bool generate_SatL_and_SatV = true);
@@ -58,7 +63,7 @@ public:
     void resize(unsigned int N);
     shared_ptr<HelmholtzEOSMixtureBackend> 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);
 
     // ***************************************************************
     // ***************************************************************

src/Backends/Incompressible/IncompressibleBackend.cpp

@@ -37,7 +37,7 @@ IncompressibleBackend::IncompressibleBackend(const std::vector<std::string> &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");
     //}

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
     /**

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,

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
 				);

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<long double> &mole_fractions){};
     void set_mass_fractions(const std::vector<long double> &mass_fractions){};
     

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);

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;
 

src/Tests/CoolProp-Tests.cpp

@@ -221,7 +221,7 @@ class TransportValidationFixture
 protected:
     long double actual, x1, x2;
     shared_ptr<CoolProp::AbstractState> 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<CoolProp::AbstractState> 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);