HS works now, for most fluids over the entire single-phase part aside for pseudo-pures which have issues

Closes https://github.com/CoolProp/CoolProp/issues/172
Closes https://github.com/CoolProp/CoolProp/issues/121

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

src/Backends/Helmholtz/FlashRoutines.cpp

@@ -136,15 +136,23 @@ 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)
+void FlashRoutines::HQ_flash(HelmholtzEOSMixtureBackend &HEOS, long double Tguess)
 {
+    SaturationSolvers::saturation_PHSU_pure_options options;
+    options.use_logdelta = false;
+    HEOS.specify_phase(iphase_twophase);
+    if (Tguess < 0){
+        options.use_guesses = true;
+        options.T = Tguess;
+        CoolProp::SaturationAncillaryFunction &rhoL = HEOS.get_components()[0]->ancillaries.rhoL;
+        CoolProp::SaturationAncillaryFunction &rhoV = HEOS.get_components()[0]->ancillaries.rhoV;
+        options.rhoL = rhoL.evaluate(Tguess);
+        options.rhoV = rhoV.evaluate(Tguess);
+    }
     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();
@@ -159,7 +167,14 @@ void FlashRoutines::QS_flash(HelmholtzEOSMixtureBackend &HEOS)
 {
     if (HEOS.is_pure_or_pseudopure){
         
-        if (std::abs(HEOS.Q()) < 1e-10){
+        if (std::abs(HEOS.smolar() - HEOS.get_state("reducing").smolar) < 0.001)
+        {
+            HEOS._p = HEOS.p_critical();
+            HEOS._T = HEOS.T_critical();
+            HEOS._rhomolar = HEOS.rhomolar_critical();
+            HEOS._phase = iphase_critical_point;
+        }
+        else 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;
@@ -1141,6 +1156,8 @@ void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, parameters ot
         HEOS.calc_pressure();
         HEOS._Q = -1;
     }
+    // Update the state for conditions where the state was guessed
+    
 }
 
 void FlashRoutines::HS_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, long double hmolar_spec, long double smolar_spec, HS_flash_singlephaseOptions &options)
@@ -1194,7 +1211,7 @@ void FlashRoutines::HS_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, long
             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);
+    while(std::abs(resid) > 1e-9);
 }
 void FlashRoutines::HS_flash_generate_TP_singlephase_guess(HelmholtzEOSMixtureBackend &HEOS, double &T, double &p)
 {
@@ -1219,47 +1236,71 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS)
         
         // Find maxima states if needed
         // Cache the maximum enthalpy saturation state;
-        //HEOS.calc_hsat_max();
+        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();
+        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
+                              
+        double first_maxima_in_saturation_entropy;
+        if (HEOS.ssat_max.exists == SsatSimpleState::SSAT_MAX_DOES_EXIST){
+            first_maxima_in_saturation_entropy = HEOS.ssat_max.smolar;
+        }
+        else{
+            first_maxima_in_saturation_entropy = tripleV.smolar;
+        }
+        
+        double h1 = HEOS.hmolar(), s1 = HEOS.smolar();
+        // Enthalpy at solid line for given entropy
         double hsolid = (tripleV.hmolar-tripleL.hmolar)/(tripleV.smolar-tripleL.smolar)*(HEOS.smolar()-tripleL.smolar) + tripleL.hmolar;
         // Part A - first check if HS is below triple line formed by connecting the triple point states
         // If so, it is solid, and not supported
-        if (HEOS.hmolar() < hsolid){
-            throw ValueError(format("Enthalpy [%g J/mol] is below solid enthalpy [%g J/mol]", HEOS.hmolar(), hsolid));
+        if (HEOS.hmolar() < hsolid-0.1){ // -0.1 is for a buffer
+            throw ValueError(format("Enthalpy [%g J/mol] is below solid enthalpy [%g J/mol] for entropy [%g J/mol/K]", HEOS.hmolar(), hsolid-0.1, HEOS.smolar()));
         }
-        /*// 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));
+        /* Now check up to the first maxima in saturated vapor entropy.
+         * For nicely behaved fluids, this means all the way up to the triple point vapor
+         * For not-nicely behaved fluids with a local maxima on the saturated vapor entropy curve, 
+         * the entropy at the local maxima in entropy
+         * 
+         * In this section, there can only be one solution for the saturation temperature, which is why the solver
+         * is divided up in this way.
+         */
+        else if (HEOS.smolar() < first_maxima_in_saturation_entropy){
+            double Q;
+            if (HEOS.smolar() < crit.smolar){ 
+                Q = 0; // Liquid part
+            }
+            else{
+                Q = 1; // Vapor part
             }
-        }*/
-        // 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){
-            
-            // Update the temporary instance with saturated liquid entropy 
-            HEOS_copy->update(QS_INPUTS, 0, HEOS.smolar());
             
+            // Update the temporary instance with saturated entropy 
+            HEOS_copy->update(QSmolar_INPUTS, Q, HEOS.smolar());
+                
             // Check if above the saturation enthalpy for given entropy
+            // If it is, the inputs are definitely single-phase.  We are done here
+            double h1 = HEOS.hmolar(), h2 = HEOS_copy->hmolar();
             if (HEOS.hmolar() > HEOS_copy->hmolar()){
                 solution = single_phase_solution;
             }
             else{
-                // C2: It is below hsatL(ssatL)
+                // C2: It is below hsat(ssat)
                 // 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())
-                {
+                /* If enthalpy is between enthalpy at maxima in h and triple
+                 * point enthalpy, search in enthalpy to find possible solutions
+                 * that yield the correct enthalpy
+                 * 
+                 * There should only be one solution since we have already bypassed the local maxima in enthalpy
+                 */
+                
+                HEOS_copy->update_HmolarQ_with_guessT(HEOS.hmolar(), 1, HEOS.hsat_max.T);
+                
+                if (HEOS.smolar() > HEOS_copy->smolar()){
                     solution = single_phase_solution;
                 }
                 else{
@@ -1268,22 +1309,6 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS)
                 }
             }
         }
-        // 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;
@@ -1311,12 +1336,22 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS)
                     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;
+                    try{
+                        // Trying again with another guessed value
+                        HEOS_copy->update(DmolarT_INPUTS, HEOS.rhomolar_critical()*1.3, HEOS.Tmax());
+                        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){
+                        // Trying again with another guessed value
+                        HEOS_copy->update(DmolarT_INPUTS, HEOS.rhomolar_critical(), 0.5*HEOS.Tmax() + 0.5*HEOS.T_critical());
+                        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:

src/Backends/Helmholtz/FlashRoutines.h

@@ -42,7 +42,8 @@ public:
     
     /// Flash for given molar enthalpy and (molar) quality
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
-    static void HQ_flash(HelmholtzEOSMixtureBackend &HEOS);
+    /// @param Tguess (optional) The guess temperature in K to start from, ignored if < 0
+    static void HQ_flash(HelmholtzEOSMixtureBackend &HEOS, long double Tguess = -1);
     
     /// Flash for mixture given temperature or pressure and (molar) quality
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used

src/Backends/Helmholtz/Fluids/Ancillaries.cpp

@@ -114,7 +114,7 @@ double SaturationAncillaryFunction::invert(double value, double min_bound, doubl
 	try{
         // Safe to expand the domain a little bit to lower temperature, absolutely cannot exceed Tmax 
         // because then you get (negative number)^(double) which is undefined.
-		return Brent(resid,min_bound,max_bound,DBL_EPSILON,1e-12,100,errstring);
+		return Brent(resid,min_bound,max_bound,DBL_EPSILON,1e-10,100,errstring);
 	}
 	catch(std::exception &e){
 		return Secant(resid,max_bound, -0.01, 1e-12, 100, errstring);

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp

@@ -658,6 +658,20 @@ void HelmholtzEOSMixtureBackend::update_DmolarT_direct(long double rhomolar, lon
     post_update();
 }
 
+void HelmholtzEOSMixtureBackend::update_HmolarQ_with_guessT(long double hmolar, long double Q, long double Tguess)
+{
+    CoolProp::input_pairs pair = CoolProp::HmolarQ_INPUTS;
+    // Set up the state
+    pre_update(pair, hmolar, Q);
+    
+    _hmolar = hmolar;
+    _Q = Q;    
+    FlashRoutines::HQ_flash(*this, Tguess);
+    
+    // Cleanup
+    post_update();
+}
+
 void HelmholtzEOSMixtureBackend::update_TP_guessrho(long double T, long double p, long double rhomolar_guess)
 {
     CoolProp::input_pairs pair = PT_INPUTS;
@@ -784,9 +798,9 @@ 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:
+        case QSmolar_INPUTS:
             _Q = value1; _smolar = value2; FlashRoutines::QS_flash(*this); break;
-        case HQ_INPUTS:
+        case HmolarQ_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()));
@@ -1112,16 +1126,27 @@ void HelmholtzEOSMixtureBackend::calc_ssat_max(void)
             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())
+    if (!ssat_max.is_valid() && ssat_max.exists != SsatSimpleState::SSAT_MAX_DOESNT_EXIST)
     {
-        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();
+        shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(get_components()));
+        Residual resid(*HEOS_copy);
+        CoolProp::SimpleState &tripleV = HEOS_copy->get_components()[0]->triple_vapor;
+        double v1 = resid.call(hsat_max.T);
+        double v2 = resid.call(tripleV.T);
+        // If there is a sign change, there is a maxima, otherwise there is no local maxima/minima
+        if (v1*v2 < 0){
+            std::string errstr;
+            Brent(resid, hsat_max.T, tripleV.T, DBL_EPSILON, 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();
+            ssat_max.exists = SsatSimpleState::SSAT_MAX_DOES_EXIST;
+        }
+        else{
+            ssat_max.exists = SsatSimpleState::SSAT_MAX_DOESNT_EXIST;
+        }
     }
 }
 void HelmholtzEOSMixtureBackend::calc_hsat_max(void)
@@ -1141,9 +1166,10 @@ void HelmholtzEOSMixtureBackend::calc_hsat_max(void)
     };
     if (!hsat_max.is_valid())
     {
-        Residualhmax residhmax(*this);
+        shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(get_components()));
+        Residualhmax residhmax(*HEOS_copy);
         std::string errstrhmax;
-        Secant(residhmax, T_critical() - 1, 0.1, 1e-8, 30, errstrhmax);
+        Brent(residhmax, T_critical() - 0.1, HEOS_copy->Ttriple() + 1, DBL_EPSILON, 1e-8, 30, errstrhmax);
         hsat_max.T = residhmax.HEOS->T();
         hsat_max.p = residhmax.HEOS->p();
         hsat_max.rhomolar = residhmax.HEOS->rhomolar();

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h

@@ -41,7 +41,8 @@ public:
     shared_ptr<ReducingFunction> Reducing;
     ExcessTerm Excess;
     PhaseEnvelopeData PhaseEnvelope;
-    SimpleState ssat_max, hsat_max;
+    SimpleState hsat_max;
+    SsatSimpleState ssat_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
@@ -77,10 +78,21 @@ public:
     void resize(unsigned int N);
     shared_ptr<HelmholtzEOSMixtureBackend> SatL, SatV; ///<
 
+    /** \brief The standard update function
+     * @param input_pair The pair of inputs that will be provided
+     * @param value1 The first input value
+     * @param value2 The second input value
+     */
     void update(CoolProp::input_pairs input_pair, double value1, double value2);
 
+    /** \brief Update with TP and a guess for rho
+     * @param T Temperature in K
+     * @param p Pressure in Pa
+     * @param rho_guess Density in mol/m^3 guessed
+     */
     void update_TP_guessrho(long double T, long double p, long double rho_guess);
     void update_DmolarT_direct(long double rhomolar, long double T);
+    void update_HmolarQ_with_guessT(long double hmolar, long double Q, long double Tguess);
 
     /** \brief Set the components of the mixture
      * 

src/Backends/Helmholtz/VLERoutines.cpp

@@ -178,6 +178,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
 
     HEOS.calc_reducing_state();
     const SimpleState & reduce = HEOS.get_reducing_state();
+    CoolProp::SimpleState crit = HEOS.get_state("reducing");
     shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL,
                                            SatV = HEOS.SatV;
 
@@ -235,27 +236,53 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
             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()));
+                shared_ptr<HelmholtzEOSMixtureBackend> HEOS_copy(new HelmholtzEOSMixtureBackend(HEOS.get_components()));
+                HEOS_copy->update(QT_INPUTS, 1, Tmin); double hTmin = HEOS_copy->hmolar();
+                HEOS_copy->update(QT_INPUTS, 1, Tmax); double hTmax = HEOS_copy->hmolar();
+                T = (Tmax-Tmin)/(hTmax-hTmin)*(HEOS.hmolar()-hTmin) + Tmin;
 			}
         }
         else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SL)
         {
-            CoolProp::SaturationAncillaryFunction &anc = HEOS.get_components()[0]->ancillaries.sL;
+            CoolPropFluid &component = *(HEOS.get_components()[0]);
+            CoolProp::SaturationAncillaryFunction &anc = component.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); 
+            // 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()))
+            {
+                T = std::max(0.99*crit.T, crit.T-0.1);
+            }
+            else{
+                long double Tmin, Tmax, Tmin_satV;
+                HEOS.calc_Tmin_sat(Tmin, Tmin_satV);
+                Tmax = HEOS.calc_Tmax_sat();
+                // Ancillary is deltas = s - hs_anchor.s
+                // First try a conventional call
+                try{
+                    T = anc.invert(specified_value - hs_anchor.smolar, Tmin, Tmax);
+                }
+                catch(std::exception &e){
+                    try{ 
+                        T = anc.invert(specified_value - hs_anchor.smolar, Tmin - 3, Tmax + 3); 
+                    }
+                    catch(std::exception &e){
+                        double vmin = anc.evaluate(Tmin);
+                        double vmax = anc.evaluate(Tmax);
+                        if (std::abs(specified_value - hs_anchor.smolar) < std::abs(vmax)){
+                            T = Tmax - 0.1;
+                        }
+                        else{
+                            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));
+                        }
+                    }
+                }
             }
-			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");
+            CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
             class Residual : public FuncWrapper1D
             {
                 public:
@@ -267,26 +294,33 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
                 }
                 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;
+                    long double resid = s_liq + component->ancillaries.sLV.evaluate(T) - s;
+                    
+                    return resid;
                 };
             };
             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);
+            // 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, DBL_EPSILON, 1e-10, 50, errstr);
             }
-            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()));
+            catch(std::exception &e){
+                long double vmax = resid.call(Tmax);
+                // If near the critical point, use a near critical guess value for T
+                if (std::abs(specified_value - hs_anchor.smolar) < std::abs(vmax)){
+                    T = std::max(0.99*crit.T, crit.T-0.1);
+                }
+                else{
+                    shared_ptr<HelmholtzEOSMixtureBackend> HEOS_copy(new HelmholtzEOSMixtureBackend(HEOS.get_components()));
+                    HEOS_copy->update(QT_INPUTS, 1, Tmin); double sTmin = HEOS_copy->smolar();
+                    HEOS_copy->update(QT_INPUTS, 1, Tmax); double sTmax = HEOS_copy->smolar();
+                    T = (Tmax-Tmin)/(sTmax-sTmin)*(HEOS.smolar()-sTmin) + Tmin;
                 }
             }
         }
@@ -294,7 +328,9 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
         {
             throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable));
         }
-		if (T > HEOS.T_critical()-1){ T -= 1; }
+        // If T from the ancillaries is above the critical temp, this will cause failure 
+        // in ancillaries for rhoV and rhoL, decrease if needed
+        T = std::min(T, static_cast<long double>(HEOS.T_critical()-0.1));
 
         // Evaluate densities from the ancillary equations
         rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T);
@@ -306,7 +342,10 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
         // be way off, and often times negative
         SatL->update(DmolarT_INPUTS, rhoL, T);
         SatV->update(DmolarT_INPUTS, rhoV, T);
-        rhoL += -(SatL->p()-SatV->p())/SatL->first_partial_deriv(iP, iDmolar, iT);
+        double rhoL_updated = rhoL -(SatL->p()-SatV->p())/SatL->first_partial_deriv(iP, iDmolar, iT);
+        
+        // Accept the update if the liquid density is greater than the vapor density
+        if (rhoL_updated > rhoV){ rhoL = rhoL_updated; }
         
         // Update the state again with the better guess for the liquid density
         SatL->update(DmolarT_INPUTS, rhoL, T);
@@ -480,16 +519,25 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
         }
 
         v = linsolve(J, negativer);
-
-        tau += options.omega*v[0];
-
-        if (options.use_logdelta){
-            deltaL = exp(log(deltaL)+options.omega*v[1]);
-            deltaV = exp(log(deltaV)+options.omega*v[2]);
-        }
-        else{
-            deltaL += options.omega*v[1];
-            deltaV += options.omega*v[2];
+        
+        // Conditions for an acceptable step are:
+        // a) tau > 1
+        // b) rhoL > rhoV or deltaL > deltaV
+        double tau0 = tau, deltaL0 = deltaL, deltaV0 = deltaV;
+        for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1)
+        {
+            tau = tau0 + omega_local*options.omega*v[0];
+            if (options.use_logdelta){
+                deltaL = exp(log(deltaL0)+omega_local*options.omega*v[1]);
+                deltaV = exp(log(deltaV0)+omega_local*options.omega*v[2]);
+            }
+            else{
+                deltaL = deltaL0 + omega_local*options.omega*v[1];
+                deltaV = deltaV0 + omega_local*options.omega*v[2];
+            }
+            if (tau > 1 && deltaL > deltaV){
+                break;
+            }
         }
 
         rhoL = deltaL*reduce.rhomolar;
@@ -808,13 +856,21 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HE
         stepL = options.omega/DELTA*( (KV-KL)*dJV-(JV-JL)*dKV);
         stepV = options.omega/DELTA*( (KV-KL)*dJL-(JV-JL)*dKL);
 
-        if (deltaL+stepL > 1 && deltaV+stepV < 1 && deltaV+stepV > 0){
-            deltaL += stepL;  deltaV += stepV;
-            rhoL = deltaL*reduce.rhomolar; rhoV = deltaV*reduce.rhomolar;
-        }
-        else{
-            throw ValueError(format("rhosatPure_Akasaka densities are crossed"));
+        long double deltaL0 = deltaL, deltaV0 = deltaV;
+        // Conditions for an acceptable step are:
+        // a) rhoL > rhoV or deltaL > deltaV
+        for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1)
+        {
+            deltaL = deltaL0 + omega_local*stepL;
+            deltaV = deltaV0 + omega_local*stepV;
+
+            if (deltaL > 1 && deltaV < 1 && deltaV > 0){
+                break;
+            }
         }
+        
+        rhoL = deltaL*reduce.rhomolar;
+        rhoV = deltaV*reduce.rhomolar;
         iter++;
         if (iter > 100){
             throw SolutionError(format("Akasaka solver did not converge after 100 iterations"));

src/DataStructures.cpp

@@ -372,8 +372,10 @@ public:
 
 input_pair_info input_pair_list[] = {
     input_pair_info(QT_INPUTS,"QT_INPUTS","Molar quality, Temperature in K"),
-    input_pair_info(QS_INPUTS,"QS_INPUTS","Molar quality, Entropy in J/mol/K"),
-    input_pair_info(HQ_INPUTS,"HQ_INPUTS","Enthalpy in J/mol, Molar quality"),
+    input_pair_info(QSmolar_INPUTS,"QS_INPUTS","Molar quality, Entropy in J/mol/K"),
+    input_pair_info(QSmass_INPUTS,"QS_INPUTS","Molar quality, Entropy in J/kg/K"),
+    input_pair_info(HmolarQ_INPUTS,"HQ_INPUTS","Enthalpy in J/mol, Molar quality"),
+    input_pair_info(HmassQ_INPUTS,"HQ_INPUTS","Enthalpy in J/kg, Molar quality"),
     input_pair_info(PQ_INPUTS,"PQ_INPUTS","Pressure in Pa, Molar quality"),
     
     input_pair_info(PT_INPUTS, "PT_INPUTS","Pressure in Pa, Temperature in K"),

src/Tests/CoolProp-Tests.cpp

@@ -1264,6 +1264,33 @@ TEST_CASE("Triple point checks", "[triple_point]")
     }
 }
 
+TEST_CASE("Test that saturation solvers solve all the way to T = Tc", "[sat_T_to_Tc]")
+{
+    std::vector<std::string> fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),',');
+    for (std::size_t i = 0; i < fluids.size(); ++i)
+    {
+        double Tc = Props1SI(fluids[i], "Tcrit");
+        std::ostringstream ss1;
+        ss1 << "Check sat_T at Tc for " << fluids[i];
+        SECTION(ss1.str(),"")
+        {
+            double pc = PropsSI("P","T",Tc,"Q",0,fluids[i]);
+            CAPTURE(pc);
+            CHECK(ValidNumber(pc));
+        }
+        for (double j = 0.1; j > 1e-10; j /= 10)
+        {
+            std::ostringstream ss2;
+            ss2 << "Check sat_T for " << fluids[i] << " for Tc - T = " << j << " K";
+            SECTION(ss2.str(),"")
+            {
+                CHECK(ValidNumber(PropsSI("D","T",Tc-j,"Q",0,fluids[i])));
+            }
+        }
+        
+    }
+}
+
 TEST_CASE("Test that reference states are correct", "[reference_states]")
 {
     std::vector<std::string> fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),',');