Add the ability to access all the critical points in C++ and python; c…

…loses CoolProp#807

include/AbstractState.h

@@ -347,7 +347,7 @@ class AbstractState {
 
     virtual void calc_viscosity_contributions(CoolPropDbl &dilute, CoolPropDbl &initial_density, CoolPropDbl &residual, CoolPropDbl &critical){ throw NotImplementedError("calc_viscosity_contributions is not implemented for this backend"); };
     virtual void calc_conductivity_contributions(CoolPropDbl &dilute, CoolPropDbl &initial_density, CoolPropDbl &residual, CoolPropDbl &critical){ throw NotImplementedError("calc_conductivity_contributions is not implemented for this backend"); };
-
+    virtual std::vector<CriticalState> calc_all_critical_points(void){ throw NotImplementedError("calc_all_critical_points is not implemented for this backend"); };
 public:
 
     AbstractState() :_fluid_type(FLUID_TYPE_UNDEFINED), _phase(iphase_unknown), _rhospline(-_HUGE), _dsplinedp(-_HUGE), _dsplinedh(-_HUGE){ clear(); }
@@ -482,6 +482,9 @@ class AbstractState {
     double rhomolar_critical(void);
     /// Return the critical molar density in kg/m^3
     double rhomass_critical(void);
+
+    /// Return the vector of critical points, including points that are unstable or correspond to negative pressure
+    std::vector<CriticalState> all_critical_points(void){ return calc_all_critical_points(); };
 
     /// Return the reducing point temperature in K
     double T_reducing(void);

include/DataStructures.h

@@ -22,6 +22,13 @@ struct SimpleState
     }
     bool is_valid(){return ValidNumber(rhomolar) && ValidNumber(T) && ValidNumber(hmolar) && ValidNumber(p);}
 };
+
+struct CriticalState : SimpleState
+{
+    bool stable;
+    CriticalState() :stable(false){ fill(_HUGE); }
+
+};
 
 /// A modified class for the state point at the maximum saturation entropy on the vapor curve
 struct SsatSimpleState : public SimpleState

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp

@@ -818,8 +818,9 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_GWP500(void)
 CoolPropDbl HelmholtzEOSMixtureBackend::calc_T_critical(void)
 {
     if (components.size() != 1){
-        std::vector<SimpleState> critpts = find_all_critical_points();
+        std::vector<CriticalState> critpts = calc_all_critical_points();
         if (critpts.size() == 1){
+            if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
             return critpts[0].T;
         }
         else{
@@ -833,8 +834,9 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_T_critical(void)
 CoolPropDbl HelmholtzEOSMixtureBackend::calc_p_critical(void)
 {
     if (components.size() != 1){
-        std::vector<SimpleState> critpts = find_all_critical_points();
+        std::vector<CriticalState> critpts = calc_all_critical_points();
         if (critpts.size() == 1){
+            if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
             return critpts[0].p;
         }
         else{
@@ -848,8 +850,9 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_p_critical(void)
 CoolPropDbl HelmholtzEOSMixtureBackend::calc_rhomolar_critical(void)
 {
     if (components.size() != 1){
-        std::vector<SimpleState> critpts = find_all_critical_points();
+        std::vector<CriticalState> critpts = calc_all_critical_points();
         if (critpts.size() == 1){
+            if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
             return critpts[0].rhomolar;
         }
         else{
@@ -3003,7 +3006,7 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv_splined(param
     }
 }
 
-void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0)
+CoolProp::CriticalState HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0)
 {
     class Resid : public FuncWrapperND{
     public:
@@ -3073,6 +3076,19 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0)
     _critical.T = T_reducing()/x[0];
     _critical.rhomolar = x[1]*rhomolar_reducing();
     _critical.p = calc_pressure_nocache(_critical.T, _critical.rhomolar);
+
+    // First (necessary but not sufficient) check of stability
+    // All eigenvalues of M* matrix must be positive
+    Eigen::MatrixXd Mstar = MixtureDerivatives::Mstar(*this, XN_INDEPENDENT);
+    Eigen::EigenSolver<Eigen::MatrixXd> es(Mstar);
+    double min_eigenvalue = es.eigenvalues().real().minCoeff();
+
+    CriticalState critical;
+    critical.stable = (min_eigenvalue > 0); // TODO: stability analysis
+    critical.T = _critical.T;
+    critical.p = _critical.p;
+    critical.rhomolar = _critical.rhomolar;
+    return critical;
 }
 
 class OneDimObjective : public FuncWrapper1DWithTwoDerivs
@@ -3109,7 +3125,7 @@ class L0CurveTracer : public FuncWrapper1DWithDeriv
 public:
     CoolProp::HelmholtzEOSMixtureBackend &HEOS;
     double delta, tau, M1_last, R_tau, R_delta;
-    std::vector<CoolProp::SimpleState> critical_points;
+    std::vector<CoolProp::CriticalState> critical_points;
     int N_critical_points;
     Eigen::MatrixXd Lstar, adjLstar, dLstardTau, d2LstardTau2, dLstardDelta;
     L0CurveTracer(HelmholtzEOSMixtureBackend &HEOS, double tau0, double delta0) : HEOS(HEOS), delta(delta0), tau(tau0), N_critical_points(0)
@@ -3188,8 +3204,7 @@ class L0CurveTracer : public FuncWrapper1DWithDeriv
             // find the critical point and store it
             if (i > 0 && M1*M1_last < 0){
                 double rhomolar = HEOS.rhomolar_reducing()*delta_new, T = HEOS.T_reducing()/tau_new;
-                HEOS.calc_critical_point(rhomolar, T);
-                CoolProp::SimpleState crit = HEOS.get_state("critical");
+                CoolProp::CriticalState crit = HEOS.calc_critical_point(rhomolar, T);
                 critical_points.push_back(crit);
                 N_critical_points++;
                 if (get_debug_level() > 0){
@@ -3206,7 +3221,7 @@ class L0CurveTracer : public FuncWrapper1DWithDeriv
 };
 
 
-std::vector<CoolProp::SimpleState> HelmholtzEOSMixtureBackend::find_all_critical_points()
+std::vector<CoolProp::CriticalState> HelmholtzEOSMixtureBackend::calc_all_critical_points()
 {
     // Store old phase
     phases old_phase = _phase;

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h

@@ -98,8 +98,8 @@ class HelmholtzEOSMixtureBackend : public AbstractState {
     CoolPropDbl calc_second_two_phase_deriv(parameters Of, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2);
     CoolPropDbl calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, CoolPropDbl x_end);
 
-    void calc_critical_point(double rho0, double T0);
-    std::vector<CoolProp::SimpleState> find_all_critical_points();
+    CriticalState calc_critical_point(double rho0, double T0);
+    std::vector<CoolProp::CriticalState> calc_all_critical_points();
 
 	/// Calculate the phase once the state is fully calculated but you aren't sure if it is liquid or gas or ...
 	void recalculate_singlephase_phase();

wrappers/Python/CoolProp/AbstractState.pxd

@@ -16,6 +16,10 @@ cdef class PyGuessesStructure:
     cpdef public double rhomolar_liq, rhomolar_vap
     cpdef public list x, y
 
+cdef class PyCriticalState:
+    cpdef public double T, p, rhomolar, hmolar, smolar
+    cpdef public bool stable    
+
 cdef class AbstractState:
     cdef cAbstractState.AbstractState *thisptr     # hold a C++ instance which we're wrapping
     cpdef update(self, constants_header.input_pairs iInput1, double Value1, double Value2)
@@ -46,6 +50,7 @@ cdef class AbstractState:
     cpdef double rhomass_critical(self) except *
     cpdef double rhomolar_critical(self) except *
     cpdef double p_critical(self) except *
+    cpdef list all_critical_points(self)
     ## Reducing point
     cpdef double T_reducing(self) except *
     cpdef double rhomolar_reducing(self) except *

wrappers/Python/CoolProp/AbstractState.pyx

@@ -5,6 +5,9 @@ cimport constants_header
 cdef class PyPhaseEnvelopeData:
     pass
 
+cdef class PyCriticalState:
+    pass
+
 cdef class PyGuessesStructure:
     def __init__(self):
         self.T = 1e999
@@ -111,6 +114,22 @@ cdef class AbstractState:
     cpdef double p_critical(self) except *:
         """ Gets the critical pressure in Pa - wrapper of c++ function :cpapi:`CoolProp::AbstractState::p_critical` """
         return self.thisptr.p_critical()
+
+    cpdef list all_critical_points(self):
+        """ Calculate all the critical points - wrapper of c++ function :cpapi:`CoolProp::AbstractState::all_critical_points` """
+        # Get all the critical points
+        cdef vector[cAbstractState.CriticalState] critpts = self.thisptr.all_critical_points()
+        cdef PyCriticalState pypt
+        cdef list collection = []
+        # Convert to python
+        for pt in critpts:
+            pypt = PyCriticalState()
+            pypt.stable = pt.stable
+            pypt.T = pt.T
+            pypt.p = pt.p
+            pypt.rhomolar = pt.rhomolar
+            collection.append(pypt)
+        return collection
     ## Reducing point
     cpdef double T_reducing(self) except *:
         """ Gets the reducing temperature in K - wrapper of c++ function :cpapi:`CoolProp::AbstractState::T_reducing` """

wrappers/Python/CoolProp/cAbstractState.pxd

@@ -10,6 +10,11 @@ cdef extern from "PhaseEnvelope.h" namespace "CoolProp":
         size_t iTsat_max, ipsat_max, icrit
         vector[double] T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap, Q
 
+cdef extern from "DataStructures.h" namespace "CoolProp":
+    cdef cppclass CriticalState:
+        double T, p, rhomolar, hmolar, smolar
+        bool stable
+
 cdef extern from "AbstractState.h" namespace "CoolProp":
 
     cdef cppclass GuessesStructure:
@@ -57,6 +62,7 @@ cdef extern from "AbstractState.h" namespace "CoolProp":
         double rhomass_critical() except +ValueError
         double rhomolar_critical() except +ValueError
         double p_critical() except +ValueError
+        vector[CriticalState] all_critical_points() except +ValueError
         ## Reducing point
         double T_reducing() except +ValueError
         double rhomolar_reducing() except +ValueError