Added first cut at PT and P&HSU solvers - much work to be done, currently only good in gas phase.

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

src/Backends/Helmholtz/FlashRoutines.cpp

@@ -1,6 +1,7 @@
 #include "VLERoutines.h"
 #include "FlashRoutines.h"
 #include "HelmholtzEOSMixtureBackend.h"
+#include "PhaseEnvelopeRoutines.h"
 
 namespace CoolProp{
 
@@ -29,8 +30,45 @@ template<class T> T dgdbeta_RachfordRice(const std::vector<T> &z, const std::vec
 
 void FlashRoutines::PT_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS)
 {
-    if (false){//HEOS.PhaseEnvelope.built){
+    if (HEOS.PhaseEnvelope.built){
         // Use the phase envelope if already constructed to determine phase boundary
+        // Determine whether you are inside (two-phase) or outside (single-phase)
+        SimpleState closest_state;
+        std::size_t i;
+        bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS, iP, HEOS._p, iT, HEOS._T, i, closest_state);
+        if (!twophase && HEOS._T > closest_state.T){
+            // Gas solution - bounded between phase envelope temperature and very high temperature
+            //
+            // Start with a guess value from SRK
+            long double rhomolar_guess = HEOS.solver_rho_Tp_SRK(HEOS._T, HEOS._p, iphase_gas);
+            
+            solver_TP_resid resid(HEOS, HEOS._T, HEOS._p);
+            std::string errstr;
+            HEOS.specify_phase(iphase_gas);
+            try{
+                // Try using Newton's method
+                long double rhomolar = Newton(resid, rhomolar_guess, 1e-10, 100, errstr);
+                // Make sure the solution is within the bounds
+                if (!is_in_closed_range(static_cast<long double>(closest_state.rhomolar), 0.0L, rhomolar)){
+                    throw ValueError("out of range");
+                }
+                HEOS.update_DmolarT_direct(rhomolar, HEOS._T);
+            }
+            catch(std::exception &e){
+                // If that fails, try a bounded solver
+                long double rhomolar = Brent(resid, closest_state.rhomolar, 1e-10, DBL_EPSILON, 1e-10, 100, errstr);
+                // Make sure the solution is within the bounds
+                if (!is_in_closed_range(static_cast<long double>(closest_state.rhomolar), 0.0L, rhomolar)){
+                    throw ValueError("out of range");
+                }
+            }
+            HEOS.unspecify_phase();
+            
+        }
+        else{
+            // Liquid solution
+            throw ValueError();
+        }
     }
     else{
         // Following the strategy of Gernert, 2014
@@ -643,7 +681,6 @@ void FlashRoutines::HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend &HE
 			// Get the value of the desired variable
 			eos = HEOS->keyed_output(other);
             pp = HEOS->p();
-            
 
 			// Difference between the two is to be driven to zero
             r = eos - value;
@@ -694,6 +731,7 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth
 {
     bool saturation_called = false;
     long double value;
+    std::cout << "PHSU" << std::endl;
     if (HEOS.imposed_phase_index != iphase_not_imposed)
     {
         // Use the phase defined by the imposed phase
@@ -701,21 +739,25 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth
     }
     else
     {
+        std::cout << "PHSU not imposed" << std::endl;
+        // Find the phase, while updating all internal variables possible
+        switch (other)
+        {
+            case iSmolar:
+                value = HEOS.smolar(); break;
+            case iHmolar:
+                value = HEOS.hmolar(); break;
+            case iUmolar:
+                value = HEOS.umolar(); break;
+            default:
+                throw ValueError(format("Input for other [%s] is invalid", get_parameter_information(other, "long").c_str()));
+        }
         if (HEOS.is_pure_or_pseudopure)
         {
 
             // Find the phase, while updating all internal variables possible
-            switch (other)
-            {
-                case iSmolar:
-                    value = HEOS.smolar(); HEOS.p_phase_determination_pure_or_pseudopure(iSmolar, value, saturation_called); break;
-                case iHmolar:
-                    value = HEOS.hmolar(); HEOS.p_phase_determination_pure_or_pseudopure(iHmolar, value, saturation_called); break;
-                case iUmolar:
-                    value = HEOS.umolar(); HEOS.p_phase_determination_pure_or_pseudopure(iUmolar, value, saturation_called); break;
-                default:
-                    throw ValueError(format("Input for other [%s] is invalid", get_parameter_information(other, "long").c_str()));
-            }
+            HEOS.p_phase_determination_pure_or_pseudopure(other, value, saturation_called);
+            
 			if (HEOS.isHomogeneousPhase())
 			{
                 // Now we use the single-phase solver to find T,rho given P,Y using a 
@@ -765,7 +807,30 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth
         }
         else
         {
-            throw NotImplementedError("HSU_P_flash does not support mixtures (yet)");
+            std::cout << format("PHSU flash for mixture\n");
+            if (HEOS.PhaseEnvelope.built){
+                // Determine whether you are inside or outside
+                SimpleState closest_state;
+                std::size_t iclosest;
+                std::cout << format("pre is inside\n");
+                bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS, iP, HEOS._p, other, value, iclosest, closest_state);
+                std::cout << format("post is inside\n");
+                
+                std::string errstr;
+                if (!twophase){
+                    PY_singlephase_flash_resid resid(HEOS, HEOS._p, other, value);
+                    // If that fails, try a bounded solver
+                    long double rhomolar = Brent(resid, closest_state.T+10, 1000, DBL_EPSILON, 1e-10, 100, errstr);
+                    HEOS.unspecify_phase();
+                }
+                else{
+                    throw ValueError("two-phase solution for Y");
+                }
+                
+            }
+            else{
+                throw ValueError("phase envelope must be built");
+            }
         }
     }    
 }

src/Backends/Helmholtz/FlashRoutines.h

@@ -13,6 +13,7 @@ state is based.
 #define FLASHROUTINES_H
 
 #include "HelmholtzEOSMixtureBackend.h"
+#include "Solvers.h"
 
 namespace CoolProp{
 
@@ -74,5 +75,87 @@ public:
     static void PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, parameters other);
 };
 
+
+/** A residual function for the rho(T,P) solver
+ */
+class solver_TP_resid : public FuncWrapper1D
+{
+public:
+    long double T, p, r, peos, rhomolar, rhor, tau, R_u, delta, dalphar_dDelta;
+    HelmholtzEOSMixtureBackend *HEOS;
+
+    solver_TP_resid(HelmholtzEOSMixtureBackend &HEOS, long double T, long double p){
+        this->HEOS = &HEOS; this->T = T; this->p = p; this->rhor = HEOS.get_reducing_state().rhomolar;
+        this->tau = HEOS.get_reducing_state().T/T; this->R_u = HEOS.gas_constant();
+    };
+    double call(double rhomolar){
+        this->rhomolar = rhomolar;
+        delta = rhomolar/rhor; // needed for derivative
+        HEOS->update_DmolarT_direct(rhomolar, T);
+        peos = HEOS->p();
+        r = (peos-p)/p;
+        return r;
+    };
+    double deriv(double rhomolar){
+        // dp/drho|T / pspecified
+        return R_u*T*(1+2*delta*HEOS->dalphar_dDelta()+pow(delta, 2)*HEOS->d2alphar_dDelta2())/p;
+    };
+};
+
+/** A residual function for the f(P, Y) solver
+ */
+class PY_singlephase_flash_resid : public FuncWrapper1D
+{
+public:
+
+    HelmholtzEOSMixtureBackend *HEOS;
+    long double p;
+    parameters other;
+    long double r, eos, value, T, rhomolar;
+    
+    int iter;
+    long double r0, r1, T1, T0, eos0, eos1, pp;
+    PY_singlephase_flash_resid(HelmholtzEOSMixtureBackend &HEOS, long double p, parameters other, long double value) : 
+            HEOS(&HEOS), p(p), other(other), value(value)
+            {
+                iter = 0;
+                // Specify the state to avoid saturation calls, but only if phase is subcritical
+                if (HEOS.phase() == iphase_liquid || HEOS.phase() == iphase_gas ){
+                    HEOS.specify_phase(HEOS.phase());
+                }
+            };
+    double call(double T){
+
+        this->T = T;
+
+        // Run the solver with T,P as inputs;
+        HEOS->update(PT_INPUTS, p, T);
+        
+        rhomolar = HEOS->rhomolar();
+        HEOS->update(DmolarT_INPUTS, rhomolar, T);
+        // Get the value of the desired variable
+        eos = HEOS->keyed_output(other);
+        pp = HEOS->p();
+
+        // Difference between the two is to be driven to zero
+        r = eos - value;
+        
+        // Store values for later use if there are errors
+        if (iter == 0){ 
+            r0 = r; T0 = T; eos0 = eos;
+        }
+        else if (iter == 1){
+            r1 = r; T1 = T; eos1 = eos; 
+        }
+        else{
+            r0 = r1; T0 = T1; eos0 = eos1;
+            r1 = r;  T1 = T; eos1 = eos;
+        }
+
+        iter++;
+        return r;
+    };
+};
+
 } /* namespace CoolProp */
 #endif /* FLASHROUTINES_H */

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp

@@ -1696,11 +1696,19 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double
     }
 
     try{
+        
         // First we try with Newton's method with analytic derivative
         double rhomolar = Newton(resid, rhomolar_guess, 1e-8, 100, errstring);
         if (!ValidNumber(rhomolar)){
             throw ValueError();
         }
+//        double rr = this->rhomolar();
+//        for (double rho = rr*0.1; rho < 1.1*rr; rho += 100){
+//            specify_phase(iphase_gas);
+//            this->update(DmolarT_INPUTS, rho, T);
+//            unspecify_phase();
+//            std::cout << format("%g %g\n", rho, this->p());
+//        }
         return rhomolar;
     }
     catch(std::exception &)
@@ -1710,6 +1718,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double
             double rhomolar = Secant(resid, rhomolar_guess, 1.1*rhomolar_guess, 1e-8, 100, errstring);
             if (!ValidNumber(rhomolar)){throw ValueError();}
             return rhomolar;
+            
         }
         catch(std::exception &)
         {

src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp

@@ -52,7 +52,9 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
     IO.x[0] = 0.6689704673;
     IO.x[1] = 0.3310295327;
     */
-    IO.rhomolar_liq *= 1.2;
+    
+    //IO.rhomolar_liq *= 1.2;
+    
     IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::P_IMPOSED;
 
     NR.call(HEOS, IO.y, IO.x, IO);
@@ -197,11 +199,12 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
         factor = std::max(factor, static_cast<long double>(1.01));
         
         // Stop if the pressure is below the starting pressure
-        if (iter > 4 && IO.p < env.p[0]){ return; }
+        if (iter > 4 && IO.p < env.p[0]){ env.built = true; std::cout << format("envelope built.\n"); return; }
         
         // Reset the failure counter
         failure_count = 0;
     }
+    
 }
 
 void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS)
@@ -282,6 +285,7 @@ void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS)
                 };
                 double call(double rhomolar_vap){
                     PhaseEnvelopeData &env = HEOS->PhaseEnvelope;
+                    IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::RHOV_IMPOSED;
                     IO.bubble_point = false;
                     IO.rhomolar_vap = rhomolar_vap;
                     IO.y = HEOS->get_mole_fractions();
@@ -329,6 +333,109 @@ void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS)
     }
 }
 
+std::vector<std::pair<std::size_t, std::size_t> > PhaseEnvelopeRoutines::find_intersections(HelmholtzEOSMixtureBackend &HEOS, parameters iInput, long double value)
+{
+    std::vector<std::pair<std::size_t, std::size_t> > intersections;
+    
+    PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
+    for (std::size_t i = 0; i < env.p.size()-1; ++i){
+        bool matched = false;
+        switch(iInput){
+            case iP:
+                if (is_in_closed_range(env.p[i], env.p[i+1], value)){ matched = true; } break;
+            case iT:
+                if (is_in_closed_range(env.T[i], env.T[i+1], value)){ matched = true; } break;
+            case iHmolar:
+                if (is_in_closed_range(env.hmolar_vap[i], env.hmolar_vap[i+1], value)){ matched = true; } break;
+            case iSmolar:
+                if (is_in_closed_range(env.smolar_vap[i], env.smolar_vap[i+1], value)){ matched = true; } break;
+            default:
+                throw ValueError(format("bad index to find_intersections"));
+        }
+        
+        if (matched){
+            std::cout << i << " " << i+1 << std::endl;
+            intersections.push_back(std::pair<std::size_t, std::size_t>(i, i+1)); 
+        }
+    }
+    return intersections;
+}
+bool PhaseEnvelopeRoutines::is_inside(HelmholtzEOSMixtureBackend &HEOS, parameters iInput1, long double value1, parameters iInput2, long double value2, std::size_t &iclosest, SimpleState &closest_state)
+{
+    PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
+    std::vector<std::pair<std::size_t, std::size_t> > intersections = find_intersections(HEOS, iInput1, value1);
+    
+    // For now, first input must be p
+    if (iInput1 != iP){throw ValueError("For now, first input must be p in is_inside");}
+    
+    // If number of intersections is 0, input is out of range, quit
+    if (intersections.size() == 0){ throw ValueError("Input is out of range; no intersections found"); }
+    
+    // If number of intersections is 1, input will be determined based on the single intersection
+    // Need to know if values increase or decrease to the right of the intersection point
+    if (intersections.size()%2 != 0){ throw ValueError("Input is weird; odd number of intersections found"); }
+    
+    // If number of intersections is even, might be a bound
+    if (intersections.size()%2 == 0){
+        if (intersections.size() != 2){throw ValueError("for now only even value accepted is 2"); }
+        std::vector<std::size_t> other_indices(4, 0);
+        std::vector<long double> *y;
+        std::vector<long double> other_values(4, 0);
+        other_indices[0] = intersections[0].first; other_indices[1] = intersections[0].second;
+        other_indices[2] = intersections[1].first; other_indices[3] = intersections[1].second;
+        
+        switch(iInput2){
+            case iT: y = &(env.T); break;
+            case iP: y = &(env.p); break;
+            case iHmolar: y = &(env.hmolar_vap); break; 
+            case iSmolar: y = &(env.smolar_vap); break;
+            default: break;
+        }
+        
+        other_values[0] = (*y)[other_indices[0]]; other_values[1] = (*y)[other_indices[1]];
+        other_values[2] = (*y)[other_indices[2]]; other_values[3] = (*y)[other_indices[3]];
+        
+        long double min_other = *(std::min_element(other_values.begin(), other_values.end()));
+        long double max_other = *(std::max_element(other_values.begin(), other_values.end()));
+        
+        std::cout << format("min: %Lg max: %Lg val: %Lg\n", min_other, max_other, value2);
+        
+        // If by using the outer bounds of the second variable, we are outside the range, 
+        // then the value is definitely not inside the phase envelope and we don't need to 
+        // do any more analysis.
+        if (!is_in_closed_range(min_other, max_other, value2)){
+            std::vector<long double> d(4, 0);
+            d[0] = std::abs(other_values[0]-value2); d[1] = std::abs(other_values[1]-value2);
+            d[2] = std::abs(other_values[2]-value2); d[3] = std::abs(other_values[3]-value2);
+            
+            // Index of minimum distance in the other_values vector
+            std::size_t idist = std::distance(d.begin(), std::min_element(d.begin(), d.end()));
+            // Index of closest point in the phase envelope
+            std::size_t iclosest = other_indices[idist];
+            
+            // Get the state for the point which is closest to the desired value - this
+            // can be used as a bounding value in the outer single-phase flash routine
+            // since you know (100%) that it is a good bound
+            closest_state.T = env.T[iclosest];
+            closest_state.p = env.p[iclosest];
+            closest_state.rhomolar = env.rhomolar_vap[iclosest];
+            closest_state.hmolar = env.hmolar_vap[iclosest];
+            closest_state.smolar = env.smolar_vap[iclosest];
+            closest_state.Q = env.Q[iclosest];
+            
+            std::cout << format("it is not inside") << std::endl;
+            return false;
+        }
+        else{
+            // Now we have to do a saturation flash call in order to determine whether or not we are inside the phase envelope or not
+            
+            // First we can interpolate using the phase envelope to get good guesses for the necessary values
+            
+            throw ValueError("For now can't be inside");
+        }
+    }
+}
+
 } /* namespace CoolProp */
 
 #endif

src/Backends/Helmholtz/PhaseEnvelopeRoutines.h

@@ -3,12 +3,49 @@
 namespace CoolProp{
 
 class PhaseEnvelopeRoutines{
-public:
+    public:
+    /** \brief Build the phase envelope
+     *
+     * @param HEOS The HelmholtzEOSMixtureBackend instance to be used
+     */
     static void build(HelmholtzEOSMixtureBackend &HEOS);
     
     /** \brief Finalize the phase envelope and calculate maxima values, critical point, etc.
+     * 
+     * @param HEOS The HelmholtzEOSMixtureBackend instance to be used
      */
     static void finalize(HelmholtzEOSMixtureBackend &HEOS);
+    
+    /** \brief Refine the phase envelope, adding points in places that are sparse
+     *
+     * @param HEOS The HelmholtzEOSMixtureBackend instance to be used
+     */
+    //static void refine(HelmholtzEOSMixtureBackend &HEOS);
+    
+    /** \brief Determine which indices bound a given value
+     * 
+     * If you provide pressure for instance, it will return each of the indices 
+     * that bound crossings in the pressure versus rhov curve.  Thus this information
+     * can be used to determine whether another input is "inside" or "outside" the phase
+     * boundary.
+     * 
+     * @param HEOS The HelmholtzEOSMixtureBackend instance to be used
+     * @param iInput The key for the variable type that is to be checked
+     * @param value The value associated with iInput
+     */
+    static std::vector<std::pair<std::size_t, std::size_t> > find_intersections(HelmholtzEOSMixtureBackend &HEOS, parameters iInput, long double value);
+    
+    /** \brief Determine whether a pair of inputs is inside or outside the phase envelope
+     * 
+     * @param HEOS The HelmholtzEOSMixtureBackend instance to be used
+     * @param iInput1 The key for the first input
+     * @param value1 The value of the first input
+     * @param iInput2 The key for the second input
+     * @param value2 The value of the second input
+     * @param iclosest The index of the phase envelope for the closest point
+     * @param closest_state A SimpleState corresponding to the closest point found on the phase envelope
+     */
+    static bool is_inside(HelmholtzEOSMixtureBackend &HEOS, parameters iInput1, long double value1, parameters iInput2, long double value2, std::size_t &iclosest, SimpleState &closest_state);
 };
     
 } /* namespace CoolProp */

src/Backends/Helmholtz/VLERoutines.cpp

@@ -972,7 +972,7 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
             throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]",IO.Nstep_max));
         }
     }
-    while(this->error_rms > 1e-9 && min_rel_change > 10*DBL_EPSILON && iter < IO.Nstep_max);
+    while(this->error_rms > 1e-9 && min_rel_change > 1000*DBL_EPSILON && iter < IO.Nstep_max);
 
     IO.Nsteps = iter;
     IO.p = p;
@@ -1002,7 +1002,6 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
     else{
         // Liquid is incipient phase, set its composition
         rSatL.set_mole_fractions(x);
-        rSatV.set_mole_fractions(y);
     }
     
     if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
@@ -1014,9 +1013,9 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
         rSatV.update_TP_guessrho(T, p, rhomolar_vap);
     }
     else{
-        throw ValueError();
+        throw ValueError("imposed variable not set for NR VLE");
     }
-
+    
     // For diagnostic purposes calculate the pressures (no derivatives are evaluated)
     long double p_liq = rSatL.p();
     long double p_vap = rSatV.p();
@@ -1029,7 +1028,6 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
     x_N_dependency_flag xN_flag = XN_DEPENDENT;
     
     if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
-        x_N_dependency_flag xN_flag = XN_DEPENDENT;
         // For the residuals F_i (equality of fugacities)
         for (std::size_t i = 0; i < N; ++i)
         {
@@ -1101,6 +1099,7 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
     error_rms = sqrt(error_rms); // Square-root (The R in RMS)
     
     // Calculate derivatives along phase boundary;
+    // Gernert thesis 3.96 and 3.97
     long double dQ_dPsat = 0, dQ_dTsat = 0;
     for (std::size_t i = 0; i < N; ++i)
     {