Added first cut at HS solver - mostly works but some failures still

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>

src/Backends/Helmholtz/FlashRoutines.cpp

@@ -136,7 +136,62 @@ void FlashRoutines::PT_flash(HelmholtzEOSMixtureBackend &HEOS)
 	HEOS._rhomolar = HEOS.solver_rho_Tp(HEOS._T, HEOS._p);
 	HEOS._Q = -1;
 }
-
+void FlashRoutines::HQ_flash(HelmholtzEOSMixtureBackend &HEOS)
+{
+    if (HEOS.is_pure_or_pseudopure){
+        if (std::abs(HEOS.Q()-1) > 1e-10){throw ValueError(format("non-unity quality not currently allowed for HQ_flash"));}
+        // Do a saturation call for given h for vapor, first with ancillaries, then with full saturation call
+        SaturationSolvers::saturation_PHSU_pure_options options;
+        options.specified_variable = SaturationSolvers::saturation_PHSU_pure_options::IMPOSED_HV;
+        options.use_logdelta = false;
+        HEOS.specify_phase(iphase_twophase);
+        SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS.hmolar(), options);
+        HEOS._p = HEOS.SatV->p();
+        HEOS._T = HEOS.SatV->T();
+        HEOS._rhomolar = HEOS.SatV->rhomolar();
+        HEOS._phase = iphase_twophase;
+    }
+    else{
+        throw NotImplementedError("QS_flash not ready for mixtures");
+    }
+}
+void FlashRoutines::QS_flash(HelmholtzEOSMixtureBackend &HEOS)
+{
+    if (HEOS.is_pure_or_pseudopure){
+        
+        if (std::abs(HEOS.Q()) < 1e-10){
+            // Do a saturation call for given s for liquid, first with ancillaries, then with full saturation call
+            SaturationSolvers::saturation_PHSU_pure_options options;
+            options.specified_variable = SaturationSolvers::saturation_PHSU_pure_options::IMPOSED_SL;
+            options.use_logdelta = false;
+            HEOS.specify_phase(iphase_twophase);
+            SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS.smolar(), options);
+            HEOS._p = HEOS.SatL->p();
+            HEOS._T = HEOS.SatL->T();
+            HEOS._rhomolar = HEOS.SatL->rhomolar();
+            HEOS._phase = iphase_twophase;
+        }
+        else if (std::abs(HEOS.Q()-1) < 1e-10)
+        {
+            // Do a saturation call for given s for vapor, first with ancillaries, then with full saturation call
+            SaturationSolvers::saturation_PHSU_pure_options options;
+            options.specified_variable = SaturationSolvers::saturation_PHSU_pure_options::IMPOSED_SV;
+            options.use_logdelta = false;
+            HEOS.specify_phase(iphase_twophase);
+            SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS.smolar(), options);
+            HEOS._p = HEOS.SatV->p();
+            HEOS._T = HEOS.SatV->T();
+            HEOS._rhomolar = HEOS.SatV->rhomolar();
+            HEOS._phase = iphase_twophase;
+        }
+        else{
+            throw ValueError(format("non-zero or 1 quality not currently allowed for QS_flash"));
+        }
+    }
+    else{
+        throw NotImplementedError("QS_flash not ready for mixtures");
+    }
+}
 void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
 {
     long double T = HEOS._T;
@@ -916,11 +971,13 @@ void FlashRoutines::HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend &HE
         // Determine why you were out of range if you can
         // 
         long double eos0 = resid.eos0, eos1 = resid.eos1;
+        std::string name = get_parameter_information(other,"short");
+        std::string units = get_parameter_information(other,"units");
         if (eos1 > eos0 && value > eos1){
-            throw ValueError(format("HSU_P_flash_singlephase_Brent could not find a solution because %s is above the maximum value of %0.12Lg", get_parameter_information(other,"short").c_str(), eos1));
+            throw ValueError(format("HSU_P_flash_singlephase_Brent could not find a solution because %s [%Lg %s] is above the maximum value of %0.12Lg %s", name.c_str(), value, units.c_str(), eos1, units.c_str()));
         }
         if (eos1 > eos0 && value < eos0){
-            throw ValueError(format("HSU_P_flash_singlephase_Brent could not find a solution because %s is below the minimum value of %0.12Lg", get_parameter_information(other,"short").c_str(), eos0));
+            throw ValueError(format("HSU_P_flash_singlephase_Brent could not find a solution because %s [%Lg %s] is below the minimum value of %0.12Lg %s", name.c_str(), value, units.c_str(), eos0, units.c_str()));
         }
         throw;
     }
@@ -1086,6 +1143,66 @@ void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, parameters ot
     }
 }
 
+void FlashRoutines::HS_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, long double hmolar_spec, long double smolar_spec, HS_flash_singlephaseOptions &options)
+{
+    int iter = 0;
+    double resid = 9e30, resid_old = 9e30;
+    CoolProp::SimpleState reducing = HEOS.get_state("reducing");
+    do{
+        // Independent variables are T0 and rhomolar0, residuals are matching h and s
+        Eigen::Vector2d r;
+        Eigen::Matrix2d J;
+        r(0) = HEOS.hmolar() - hmolar_spec;
+        r(1) = HEOS.smolar() - smolar_spec;
+        J(0,0) = HEOS.first_partial_deriv(iHmolar, iTau, iDelta);
+        J(0,1) = HEOS.first_partial_deriv(iHmolar, iDelta, iTau);
+        J(1,0) = HEOS.first_partial_deriv(iSmolar, iTau, iDelta);
+        J(1,1) = HEOS.first_partial_deriv(iSmolar, iDelta, iTau);
+        // Step in v obtained from Jv = -r
+        Eigen::Vector2d v = J.colPivHouseholderQr().solve(-r);
+        bool good_solution = false;
+        double tau0 = HEOS.tau(), delta0 = HEOS.delta();
+        // Calculate the old residual after the last step
+        resid_old = sqrt(POW2(HEOS.hmolar() - hmolar_spec) + POW2(HEOS.smolar() - smolar_spec));
+        for (double frac = 1.0; frac > 0.001; frac /= 2)
+        {
+            try{
+                // Calculate new values
+                double tau_new = tau0 + options.omega*frac*v(0);
+                double delta_new = delta0 + options.omega*frac*v(1);
+                double T_new = reducing.T/tau_new;
+                double rhomolar_new = delta_new*reducing.rhomolar;
+                // Update state with step
+                HEOS.update(DmolarT_INPUTS, rhomolar_new, T_new);
+                resid = sqrt(POW2(HEOS.hmolar() - hmolar_spec) + POW2(HEOS.smolar() - smolar_spec));
+                if (resid > resid_old){
+                    throw ValueError(format("residual not decreasing; frac: %g, resid: %g, resid_old: %g", frac, resid, resid_old));
+                }
+                good_solution = true;
+                break;
+            }
+            catch(std::exception &e){
+                HEOS.clear();
+                continue;
+            }            
+        }
+        if (!good_solution){
+            throw ValueError(format("Not able to get a solution" ));
+        }
+        iter++;
+        if (iter > 50){
+            throw ValueError(format("HS_flash_singlephase took too many iterations; residual is %g; prior was %g", resid, resid_old));
+        }
+    }
+    while(std::abs(resid) > 1e-10);
+}
+void FlashRoutines::HS_flash_generate_TP_singlephase_guess(HelmholtzEOSMixtureBackend &HEOS, double &T, double &p)
+{
+    // Randomly obtain a starting value that is single-phase
+    double logp = ((double)rand()/(double)RAND_MAX)*(log(HEOS.pmax())-log(HEOS.p_triple()))+log(HEOS.p_triple());
+    T = ((double)rand()/(double)RAND_MAX)*(HEOS.Tmax()-HEOS.Ttriple())+HEOS.Ttriple();
+    p = exp(logp);
+}
 void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS)
 {
     if (HEOS.imposed_phase_index != iphase_not_imposed)
@@ -1095,7 +1212,19 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS)
     }
     else
     {
-        CoolProp::SimpleState &crit = HEOS.components[0]->pEOS->reduce;
+        enum solution_type_enum{not_specified = 0, single_phase_solution, two_phase_solution};
+        solution_type_enum solution;
+        
+        shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(HEOS.components));
+        
+        // Find maxima states if needed
+        // Cache the maximum enthalpy saturation state;
+        //HEOS.calc_hsat_max();
+        // For weird fluids like the siloxanes, there can also be a maximum 
+        // entropy along the vapor saturation line. Try to find it if it has one
+        // HEOS.calc_ssat_max();
+        
+        CoolProp::SimpleState crit = HEOS.get_state("reducing");
         CoolProp::SimpleState &tripleL = HEOS.components[0]->triple_liquid,
                               &tripleV = HEOS.components[0]->triple_vapor;
         // Enthalpy at solid line
@@ -1105,57 +1234,99 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS)
         if (HEOS.hmolar() < hsolid){
             throw ValueError(format("Enthalpy [%g J/mol] is below solid enthalpy [%g J/mol]", HEOS.hmolar(), hsolid));
         }
-        // Part B - Check lower limit
+        /*// Part B - Check lower limit
         else if (HEOS.smolar() < tripleL.smolar){
             // If fluid is other than water (which can have solutions below tripleL), cannot have any solution, fail
             if (upper(HEOS.name()) != "Water"){
                 throw ValueError(format("Entropy [%g J/mol/K] is below triple point liquid entropy [%g J/mol/K]", HEOS.smolar(), tripleL.smolar));
             }
-        }
-        // Part C - Check higher limit
-        else if (HEOS.smolar() > tripleV.smolar){
-            throw ValueError(format("Entropy [%g J/mol/K] is above triple point vapor entropy [%g J/mol/K]", HEOS.smolar(), tripleV.smolar));
-        }
-        // Part D - if s < sc, a few options are possible.  It could be two-phase, or liquid (more likely), or gas (less likely)
+        }*/
+        // Part C - if s < sc, a few options are possible.  It could be two-phase, or liquid (more likely), or gas (less likely)
         else if (HEOS.smolar() < crit.smolar){
-            // Do a saturation call for given s for liquid, first with ancillaries, then with full saturation call
-            // Check if above the saturation enthalpy each time
-            // D1: It is above hsat
-            // It is liquid, supercritical liquid, or out of range
-            // Come up with some guess values, maybe start at the saturation curve with under-relaxation?
-            // D2: It is below hsat
-            // Either two-phase, or vapor (for funky saturation curves like the siloxanes)
-            // Do a saturation_h call for given h on the vapor line to determine whether two-phase or vapor
-            // D2a: It is below ssat
-            // two-phase
-            // D2b: It is above ssat
-            // vapor
-            throw ValueError("partD");
+            
+            // Update the temporary instance with saturated liquid entropy 
+            HEOS_copy->update(QS_INPUTS, 0, HEOS.smolar());
+            
+            // Check if above the saturation enthalpy for given entropy
+            if (HEOS.hmolar() > HEOS_copy->hmolar()){
+                solution = single_phase_solution;
+            }
+            else{
+                // C2: It is below hsatL(ssatL)
+                // Either two-phase, or vapor (for funky saturation curves like the siloxanes)
+                // Do a saturation_h call for given h on the vapor line to determine whether two-phase or vapor
+                // Update the temporary instance with saturated vapor enthalpy
+                HEOS_copy->update(HQ_INPUTS, HEOS.hmolar(), 1);
+                
+                if (HEOS.smolar() > HEOS_copy->smolar())
+                {
+                    solution = single_phase_solution;
+                }
+                else{
+                    // C2a: It is below ssatV(hsatV) --> two-phase
+                    solution = two_phase_solution;
+                }
+            }
+        }
+        // Part C - if tripleV.s > s > sc
+        else if (HEOS.smolar() > crit.smolar && HEOS.smolar() < tripleV.smolar){
+            // Do a saturation_s call for given s on the vapor line to determine whether two-phase or vapor
+            // Update the temporary instance with saturated vapor entropy
+            HEOS_copy->update(QS_INPUTS, 1, HEOS.smolar());
+            
+            double h1 = HEOS.hmolar(), h2 = HEOS_copy->hmolar();
+            if (HEOS.hmolar() > HEOS_copy->hmolar()){
+                // D2b: It is above hsatV(ssatV) --> gas
+                solution = single_phase_solution;
+            }
+            else{
+                // C2a: It is below ssatV(hsatV) --> two-phase
+                solution = two_phase_solution;
+            }
+        }
+        // Part D - Check higher limit
+        else if (HEOS.smolar() > tripleV.smolar){
+            solution = single_phase_solution;
+            HEOS_copy->update(PT_INPUTS, HEOS_copy->p_triple(), 0.5*HEOS_copy->Tmin() + 0.5*HEOS_copy->Tmax());
         }
         // Part E - HEOS.smolar() > crit.hmolar > tripleV.smolar
         else{
-            throw ValueError("partE");
+            // Now branch depending on the saturated vapor curve
+            // If maximum 
+            throw ValueError(format("partE HEOS.smolar() = %g tripleV.smolar = %g", HEOS.smolar(), tripleV.smolar));
+        }
+        
+        switch (solution){
+            case single_phase_solution:
+            {
+                // Fixing it to be gas is probably sufficient
+                HEOS_copy->specify_phase(iphase_gas);
+                HS_flash_singlephaseOptions options;
+                options.omega = 1.0;
+                try{
+                    // Do the flash calcs starting from the guess value
+                    HS_flash_singlephase(*HEOS_copy, HEOS.hmolar(), HEOS.smolar(), options);
+                    // Copy the results
+                    HEOS.update(DmolarT_INPUTS, HEOS_copy->rhomolar(), HEOS_copy->T());
+                    break;
+                }
+                catch(std::exception &e){
+                    HEOS_copy->update(DmolarT_INPUTS, HEOS.rhomolar_critical()*1.3, HEOS.Tmax());
+                    // Do the flash calcs starting from the guess value
+                    HS_flash_singlephase(*HEOS_copy, HEOS.hmolar(), HEOS.smolar(), options);
+                    // Copy the results
+                    HEOS.update(DmolarT_INPUTS, HEOS_copy->rhomolar(), HEOS_copy->T());
+                    break;
+                }
+            }
+            case two_phase_solution:
+            {
+                throw ValueError("HS two-phase not yet supported.");
+            }
+            default:
+                throw ValueError("solution not set");
         }
     }
-
-    /*if (HEOS.isHomogeneousPhase() && !ValidNumber(HEOS._p))
-    {
-        switch (other)
-        {
-            case iDmolar:
-                break;
-            case iHmolar:
-                HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._hmolar, iHmolar); break;
-            case iSmolar:
-                HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._smolar, iSmolar); break;
-            case iUmolar:
-                HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._umolar, iUmolar); break;
-            default:
-                break;
-        }
-        HEOS.calc_pressure();
-        HEOS._Q = -1;
-    }*/
 }
 
 } /* namespace CoolProp */

src/Backends/Helmholtz/FlashRoutines.h

@@ -36,6 +36,14 @@ public:
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
     static void QT_flash(HelmholtzEOSMixtureBackend &HEOS);
     
+    /// Flash for given molar entropy and (molar) quality
+    /// @param HEOS The HelmholtzEOSMixtureBackend to be used
+    static void QS_flash(HelmholtzEOSMixtureBackend &HEOS);
+    
+    /// Flash for given molar enthalpy and (molar) quality
+    /// @param HEOS The HelmholtzEOSMixtureBackend to be used
+    static void HQ_flash(HelmholtzEOSMixtureBackend &HEOS);
+    
     /// Flash for mixture given temperature or pressure and (molar) quality
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
     /// @param other The parameter that is imposed, either iT or iP
@@ -83,6 +91,19 @@ public:
     /// A flash routine for (H,S)
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
     static void HS_flash(HelmholtzEOSMixtureBackend &HEOS);
+    
+    /// Randomly generate a single phase set of inputs for T and p - searches entire single-phase region
+    /// @param HEOS The HelmholtzEOSMixtureBackend to be used
+    /// @param T The temperature in K
+    /// @param p The pressure in Pa
+    static void HS_flash_generate_TP_singlephase_guess(HelmholtzEOSMixtureBackend &HEOS, double &T, double &p);
+    
+    struct HS_flash_singlephaseOptions
+    {
+        double omega;
+        HS_flash_singlephaseOptions(){omega = 1.0;}
+    };
+    static void HS_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, long double hmolar_spec, long double smolar_spec, HS_flash_singlephaseOptions &options);
 };
 
 

src/Backends/Helmholtz/Fluids/FluidLibrary.h

@@ -313,6 +313,8 @@ protected:
         EOS.reduce.T = cpjson::get_double(reducing_state,"T");
         EOS.reduce.rhomolar = cpjson::get_double(reducing_state,"rhomolar");
         EOS.reduce.p = cpjson::get_double(reducing_state,"p");
+        EOS.reduce.hmolar = cpjson::get_double(reducing_state,"hmolar");
+        EOS.reduce.smolar = cpjson::get_double(reducing_state,"smolar");
 
         EOS.sat_min_liquid.T = cpjson::get_double(satminL_state, "T");
         EOS.sat_min_liquid.p = cpjson::get_double(satminL_state, "p");
@@ -947,7 +949,7 @@ protected:
         rapidjson::Value &crit = states["critical"];
         fluid.crit.T = cpjson::get_double(crit, "T");
         fluid.crit.p = cpjson::get_double(crit, "p");
-        fluid.crit.rhomolar = cpjson::get_double(crit, "rhomolar");
+        fluid.crit.rhomolar = cpjson::get_double(crit, "rhomolar");        
 
         if (!states.HasMember("triple_liquid")){ throw ValueError(format("fluid[\"STATES\"] [%s] does not have \"triple_liquid\" member",fluid.name.c_str())); }
         rapidjson::Value &triple_liquid = states["triple_liquid"];

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp

@@ -784,6 +784,10 @@ void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double
             _Q = value1; _T = value2; FlashRoutines::QT_flash(*this); break;
         case PQ_INPUTS:
             _p = value1; _Q = value2; FlashRoutines::PQ_flash(*this); break;
+        case QS_INPUTS:
+            _Q = value1; _smolar = value2; FlashRoutines::QS_flash(*this); break;
+        case HQ_INPUTS:
+            _hmolar = value1; _Q = value2; FlashRoutines::HQ_flash(*this); break;
         default:
             throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
     }
@@ -796,7 +800,8 @@ 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");}
+    if (!ValidNumber(_p)){ 
+        throw ValueError("p is not a valid number");}
     //if (_T < 0){ throw ValueError("T is less than zero");}
     if (!ValidNumber(_T)){ throw ValueError("T is not a valid number");}
     if (_rhomolar < 0){ throw ValueError("rhomolar is less than zero");}
@@ -1092,6 +1097,60 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
         throw ValueError(format("The pressure [%g Pa] cannot be used in p_phase_determination",_p));
     }
 }
+void HelmholtzEOSMixtureBackend::calc_ssat_max(void)
+{
+    class Residual : public FuncWrapper1D
+    {
+        public:
+        HelmholtzEOSMixtureBackend *HEOS;
+        Residual(HelmholtzEOSMixtureBackend &HEOS): HEOS(&HEOS){};
+        double call(double T){
+            HEOS->update(QT_INPUTS, 1, T);
+            // dTdp_along_sat
+            double dTdp_along_sat = HEOS->T()*(1/HEOS->SatV->rhomolar()-1/HEOS->SatL->rhomolar())/(HEOS->SatV->hmolar()-HEOS->SatL->hmolar());
+            // dsdT_along_sat;
+            return HEOS->SatV->first_partial_deriv(iSmolar,iT,iP)+HEOS->SatV->first_partial_deriv(iSmolar,iP,iT)/dTdp_along_sat;
+        }
+    };
+    if (!ssat_max.is_valid())
+    {
+        Residual resid(*this);
+        std::string errstr;
+        Secant(resid, T_critical() - 1, 0.1, 1e-8, 30, errstr);
+        ssat_max.T = resid.HEOS->T();
+        ssat_max.p = resid.HEOS->p();
+        ssat_max.rhomolar = resid.HEOS->rhomolar();
+        ssat_max.hmolar = resid.HEOS->hmolar();
+        ssat_max.smolar = resid.HEOS->smolar();
+    }
+}
+void HelmholtzEOSMixtureBackend::calc_hsat_max(void)
+{
+    class Residualhmax : public FuncWrapper1D
+    {
+        public:
+        HelmholtzEOSMixtureBackend *HEOS;
+        Residualhmax(HelmholtzEOSMixtureBackend &HEOS): HEOS(&HEOS){};
+        double call(double T){
+            HEOS->update(QT_INPUTS, 1, T);
+            // dTdp_along_sat
+            double dTdp_along_sat = HEOS->T()*(1/HEOS->SatV->rhomolar()-1/HEOS->SatL->rhomolar())/(HEOS->SatV->hmolar()-HEOS->SatL->hmolar());
+            // dhdT_along_sat;
+            return HEOS->SatV->first_partial_deriv(iHmolar,iT,iP)+HEOS->SatV->first_partial_deriv(iHmolar,iP,iT)/dTdp_along_sat;
+        }
+    };
+    if (!hsat_max.is_valid())
+    {
+        Residualhmax residhmax(*this);
+        std::string errstrhmax;
+        Secant(residhmax, T_critical() - 1, 0.1, 1e-8, 30, errstrhmax);
+        hsat_max.T = residhmax.HEOS->T();
+        hsat_max.p = residhmax.HEOS->p();
+        hsat_max.rhomolar = residhmax.HEOS->rhomolar();
+        hsat_max.hmolar = residhmax.HEOS->hmolar();
+        hsat_max.smolar = residhmax.HEOS->smolar();
+    }
+}
 void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int other, long double value)
 {
     if (!ValidNumber(value)){

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h

@@ -41,12 +41,13 @@ public:
     shared_ptr<ReducingFunction> Reducing;
     ExcessTerm Excess;
     PhaseEnvelopeData PhaseEnvelope;
+    SimpleState ssat_max, hsat_max;
 
     friend class FlashRoutines; // Allows the static methods in the FlashRoutines class to have access to all the protected members and methods of this class
     friend class TransportRoutines; // Allows the static methods in the TransportRoutines class to have access to all the protected members and methods of this class
     friend class MixtureDerivatives; // Allows the static methods in the MixtureDerivatives class to have access to all the protected members and methods of this class
     friend class PhaseEnvelopeRoutines; // Allows the static methods in the PhaseEnvelopeRoutines class to have access to all the protected members and methods of this class
-    friend class MixtureParameters; //  Allows the static methods in the MixtureParameters class to have access to all the protected members and methods of this class
+    friend class MixtureParameters; // Allows the static methods in the MixtureParameters class to have access to all the protected members and methods of this class
 
     // Helmholtz EOS backend uses mole fractions
     bool using_mole_fractions(){return true;}
@@ -59,6 +60,8 @@ public:
     void calc_specify_phase(phases phase){ specify_phase(phase); }
     void calc_unspecify_phase(){ unspecify_phase(); }
     long double calc_saturation_ancillary(parameters param, int Q, parameters given, double value);
+    void calc_ssat_max(void);
+    void calc_hsat_max(void);
     
     const CoolProp::SimpleState &calc_state(const std::string &state);
 

src/Backends/Helmholtz/VLERoutines.cpp

@@ -181,7 +181,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
     shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL,
                                            SatV = HEOS.SatV;
 
-    long double T, rhoL, rhoV, pL, pV;
+    long double T, rhoL, rhoV, pL, pV, hL, sL, hV, sV;
     long double deltaL=0, deltaV=0, tau=0, error;
     int iter=0, specified_parameter;
 
@@ -200,6 +200,96 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
 				throw ValueError("Unable to invert ancillary equation");
 			}
         }
+        else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_HL)
+        {
+            CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
+            // Ancillary is deltah = h - hs_anchor.h
+            try{ T = HEOS.get_components()[0]->ancillaries.hL.invert(specified_value - hs_anchor.hmolar); }
+			catch(std::exception &e){
+				throw ValueError("Unable to invert ancillary equation for hL");
+			}
+        }
+        else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_HV)
+        {
+            class Residual : public FuncWrapper1D
+            {
+                public:
+                CoolPropFluid *component;
+                double h;
+                Residual(CoolPropFluid &component, double h){
+                    this->component = &component;
+                    this->h = h;
+                }
+                double call(double T){
+                    long double h_liq = component->ancillaries.hL.evaluate(T) + component->pEOS->hs_anchor.hmolar;
+                    return h_liq + component->ancillaries.hLV.evaluate(T) - h;
+                };
+            };
+            Residual resid(*(HEOS.get_components()[0]), HEOS.hmolar());
+            
+            // Ancillary is deltah = h - hs_anchor.h
+            std::string errstr;
+            long double Tmin_satL, Tmin_satV;
+            HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
+            double Tmin = Tmin_satL;
+            double Tmax = HEOS.calc_Tmax_sat();
+            try{ T = Brent(resid, Tmin-3, Tmax + 1, DBL_EPSILON, 1e-10, 50, errstr); }
+			catch(std::exception &e){
+				throw ValueError(format("Unable to invert ancillary equation for hV: %s",e.what()));
+			}
+        }
+        else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SL)
+        {
+            CoolProp::SaturationAncillaryFunction &anc = HEOS.get_components()[0]->ancillaries.sL;
+            CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
+            // Ancillary is deltas = s - hs_anchor.s
+            try{ 
+                T = anc.invert(specified_value - hs_anchor.smolar); 
+            }
+			catch(std::exception &e){
+                double vmin = anc.evaluate(anc.get_Tmin());
+                double vmax = anc.evaluate(anc.get_Tmax());
+				throw ValueError(format("Unable to invert ancillary equation for sL for value %Lg with Tminval %g and Tmaxval %g ", specified_value - hs_anchor.smolar, vmin, vmax));
+			}
+        }
+        else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SV)
+        {
+            CoolPropFluid &component = *(HEOS.get_components()[0]);
+            CoolProp::SimpleState crit = HEOS.get_state("reducing");
+            class Residual : public FuncWrapper1D
+            {
+                public:
+                CoolPropFluid *component;
+                double s;
+                Residual(CoolPropFluid &component, double s){
+                    this->component = &component;
+                    this->s = s;
+                }
+                double call(double T){
+                    long double s_liq = component->ancillaries.sL.evaluate(T) + component->pEOS->hs_anchor.smolar;
+                    return s_liq + component->ancillaries.sLV.evaluate(T) - s;
+                };
+            };
+            Residual resid(component, HEOS.smolar());
+            
+            // If near the critical point, use a near critical guess value for T
+            if (std::abs(HEOS.smolar() - crit.smolar) < std::abs(component.ancillaries.sL.get_max_abs_error() + component.ancillaries.sLV.get_max_abs_error()))
+            {
+                T = std::max(0.99*crit.T, crit.T-0.1);
+            }
+            else{
+                // Ancillary is deltas = s - hs_anchor.s
+                std::string errstr;
+                long double Tmin_satL, Tmin_satV;
+                HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
+                double Tmin = Tmin_satL;
+                double Tmax = HEOS.calc_Tmax_sat();
+                try{ T = Brent(resid, Tmin-3, Tmax+1, DBL_EPSILON, 1e-10, 50, errstr); }
+                catch(std::exception &e){
+                    throw ValueError(format("Unable to invert ancillary equation for sV: %s",e.what()));
+                }
+            }
+        }
         else
         {
             throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable));
@@ -237,7 +327,11 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
         }*/
 
         pL = SatL->p();
+        hL = SatL->hmolar();
+        sL = SatL->smolar();
         pV = SatV->p();
+        hV = SatV->hmolar();
+        sV = SatV->smolar();
 
         // These derivatives are needed for both cases
         long double alpharL = SatL->alphar();
@@ -251,17 +345,29 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
         long double d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
         long double d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
 
-        // -r_1
+        // -r_1 (equate the pressures)
         negativer[0] = -(deltaV*(1+deltaV*dalphar_ddeltaV)-deltaL*(1+deltaL*dalphar_ddeltaL));
-        // -r_2
+        // -r_2 (equate the gibbs energy)
         negativer[1] = -(deltaV*dalphar_ddeltaV+alpharV+log(deltaV)-deltaL*dalphar_ddeltaL-alpharL-log(deltaL));
         switch (options.specified_variable){
             case saturation_PHSU_pure_options::IMPOSED_PL:
-                // -r_3
+                // -r_3 (equate calculated pressure and specified liquid pressure)
                 negativer[2] = -(pL/specified_value - 1); break;
             case saturation_PHSU_pure_options::IMPOSED_PV:
-                // -r_3
+                // -r_3 (equate calculated pressure and specified vapor pressure)
                 negativer[2] = -(pV/specified_value - 1); break;
+            case saturation_PHSU_pure_options::IMPOSED_HL:
+                // -r_3 (equate calculated liquid enthalpy and specified liquid enthalpy)
+                negativer[2] = -(hL - specified_value); break;
+            case saturation_PHSU_pure_options::IMPOSED_HV:
+                // -r_3 (equate calculated vapor enthalpy and specified vapor enthalpy)
+                negativer[2] = -(hV - specified_value); break;
+            case saturation_PHSU_pure_options::IMPOSED_SL:
+                // -r_3 (equate calculated liquid entropy and specified liquid entropy)
+                negativer[2] = -(sL - specified_value); break;
+            case saturation_PHSU_pure_options::IMPOSED_SV:
+                // -r_3 (equate calculated vapor entropy and specified vapor entropy)
+                negativer[2] = -(sV - specified_value); break;
             default:
                 throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable));
         }
@@ -329,6 +435,46 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
                 }
 				specified_parameter = CoolProp::iP;
                 break;
+            case saturation_PHSU_pure_options::IMPOSED_HL:
+                // dr_3/dtau
+                J[2][0] = SatL->first_partial_deriv(iHmolar,iTau,iDelta);
+                // dr_3/ddelta'
+                J[2][1] = SatL->first_partial_deriv(iHmolar,iDelta,iTau);
+                if (options.use_logdelta){ J[2][1]*=deltaL;}
+                // dr_3/ddelta''
+                J[2][2] = 0; //(liquid enthalpy not a function of vapor density)
+				specified_parameter = CoolProp::iHmolar;
+                break;
+            case saturation_PHSU_pure_options::IMPOSED_HV:
+                // dr_3/dtau
+                J[2][0] = SatV->first_partial_deriv(iHmolar,iTau,iDelta);
+                // dr_3/ddelta'
+                J[2][1] = 0; //(vapor enthalpy not a function of liquid density)
+                // dr_3/ddelta''
+                J[2][2] = SatV->first_partial_deriv(iHmolar,iDelta,iTau);
+                if (options.use_logdelta){ J[2][2]*=deltaV;}
+				specified_parameter = CoolProp::iHmolar;
+                break;
+            case saturation_PHSU_pure_options::IMPOSED_SL:
+                // dr_3/dtau
+                J[2][0] = SatL->first_partial_deriv(iSmolar,iTau,iDelta);
+                // dr_3/ddelta'
+                J[2][1] = SatL->first_partial_deriv(iSmolar,iDelta,iTau);
+                if (options.use_logdelta){ J[2][1] *= deltaL; }
+                // dr_3/ddelta''
+                J[2][2] = 0; //(liquid entropy not a function of vapor density)
+				specified_parameter = CoolProp::iSmolar;
+                break;
+            case saturation_PHSU_pure_options::IMPOSED_SV:
+                // dr_3/dtau
+                J[2][0] = SatV->first_partial_deriv(iSmolar,iTau,iDelta);
+                // dr_3/ddelta'
+                J[2][1] = 0; //(vapor enthalpy not a function of liquid density)
+                // dr_3/ddelta''
+                J[2][2] = SatV->first_partial_deriv(iSmolar,iDelta,iTau);
+                if (options.use_logdelta){ J[2][2]*=deltaV;}
+				specified_parameter = CoolProp::iSmolar;
+                break;
             default:
                 throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable));
         }
@@ -359,12 +505,12 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
         {
             throw SolutionError(format("saturation_PHSU_pure solver T < 0"));
         }
-        if (iter > 25){
+        if (iter > 50){
             // Set values back into the options structure for use in next solver
             options.rhoL = rhoL; options.rhoV = rhoV; options.T = T;
             // Error out
 			std::string info = get_parameter_information(specified_parameter, "short");
-            throw SolutionError(format("saturation_PHSU_pure solver did not converge after 25 iterations for %s=%Lg current error is %Lg", info.c_str(), specified_value, error));
+            throw SolutionError(format("saturation_PHSU_pure solver did not converge after 50 iterations for %s=%Lg current error is %Lg", info.c_str(), specified_value, error));
         }
     }
     while (error > 1e-9);

src/Tests/CoolProp-Tests.cpp

@@ -1333,6 +1333,7 @@ TEST_CASE("Test that reference states are correct", "[reference_states]")
                 double EOS_smolar = HEOS->smolar();
                 CHECK( std::abs(EOS_hmolar - reducing.hmolar) < 1e-3);
                 CHECK( std::abs(EOS_smolar - reducing.smolar) < 1e-3);
+                CHECK(ValidNumber(reducing.hmolar));
                 // Then set the reference state back to the default
                 set_reference_stateS(fluids[i],"RESET");
             }   
@@ -1351,7 +1352,7 @@ TEST_CASE("Test that HS solver works for a few fluids", "[HS_solver]")
         {
             double Tmin = HEOS->Ttriple();
             double Tmax = HEOS->Tmax();
-            for (double T = Tmin; T < Tmax; T += 100)
+            for (double T = Tmin + 1; T < Tmax-1; T += 10)
             {
                 std::ostringstream ss;
                 ss << "Check HS for " << fluids[i] << " for T=" << T << ", p=" << p;
@@ -1367,6 +1368,9 @@ TEST_CASE("Test that HS solver works for a few fluids", "[HS_solver]")
                         CAPTURE(HEOS->hmolar());
                         CAPTURE(HEOS->smolar());
                         CHECK_NOTHROW(HEOS->update(HmolarSmolar_INPUTS, HEOS->hmolar(), HEOS->smolar()));
+                        double Terr = HEOS->T()- T;
+                        CAPTURE(Terr);
+                        CHECK(std::abs(Terr) < 1e-6);
                     }
                 }
             }