enhancements for setting method and coolprop free energy

thermosteam/_chemical.py

@@ -107,6 +107,25 @@ def Z(self, T, P):
     def __repr__(self):
         return f"{type(self).__name__}(self.V)"
 
+try:
+    from CoolProp.CoolProp import PropsSI
+except:
+    pass
+else:
+    class CoolPropFreeEnergy:
+        __slots__ = ('_ref', 'cache', 'name', 'CAS', 'MW')
+
+        def __init__(self, CAS, name, MW):
+            self.CAS = CAS
+            self.name = name
+            self.cache = (None, None)
+            self._ref = PropsSI(self.name, 'T', 298.15, 'P', 101325., self.CAS) / 1000. * MW
+            self.MW = MW
+
+        def __call__(self, T, P):
+            return PropsSI(self.name, 'T', T, 'P', P, self.CAS) * self.MW / 1000. - self._ref
+
+
 # %% Filling missing properties
 
 def get_chemical_data(chemical):
@@ -487,7 +506,7 @@ def __new__(cls, ID, cache=None, *, search_ID=None,
                 default=False, phase=None, search_db=True, 
                 V=None, Cn=None, mu=None, Cp=None, rho=None,
                 sigma=None, kappa=None, epsilon=None, Psat=None,
-                Hvap=None, **data):
+                Hvap=None, method=None, **data):
         chemical_cache = cls.chemical_cache
         if (cache or cls.cache) and ID in chemical_cache:
             if any([search_ID, eos, phase_ref, CAS, default, phase, 
@@ -537,8 +556,49 @@ def __new__(cls, ID, cache=None, *, search_ID=None,
                 for i in chemical_cache:
                     del chemical_cache[i]
                     break
+        if method: self.set_method(method)
         return self
 
+    def set_method(self, method):
+        for name in _model_handles:
+            model_handle = getattr(self, name)
+            if method not in model_handle.all_methods: continue
+            model_handle.method = method
+        for name in _phase_handles:
+            obj = getattr(self, name)
+            if method in model_handle.all_methods: 
+                if isinstance(obj, PhaseHandle):
+                    for phase, model_handle in obj:
+                        model_handle.method = method
+
+                else:
+                    obj.method = method
+            if hasattr(model_handle, 'all_methods_P') and method in model_handle.all_methods_P: 
+                if isinstance(obj, PhaseHandle):
+                    for phase, model_handle in obj:
+                        model_handle.method_P = method
+                else:
+                    obj.method_P = method
+
+    def use_coolprop_free_energies(self):
+        setattr = object.__setattr__
+        for name in ('_H', '_S'):
+            obj = getattr(self, name)
+            # Phase determined by coolprop
+            cpfe = CoolPropFreeEnergy(self.CAS, name.strip('_'), self.MW)
+            if isinstance(obj, PhaseHandle):
+                for phase, model_handle in obj: setattr(obj, phase, cpfe)
+            else:
+                setattr(obj, phase, cpfe)
+        # Excess energies are always included
+        f = lambda T, P: 0.
+        for name in ('_H_excess', '_S_excess'):
+            obj = getattr(self, name)
+            if isinstance(obj, PhaseHandle):
+                for phase, model_handle in obj: setattr(obj, phase, f)
+            else:
+                setattr(obj, phase, f)
+
     @classmethod
     def new(cls, ID, CAS, eos=None, phase_ref=None, phase=None, **data):
         """Create a new chemical from data without searching through

thermosteam/base/phase_handle.py

@@ -112,7 +112,7 @@ class PhaseTPHandle(PhaseHandle):
     __slots__ = ()
     force_gas_critical_phase = False
 
-    def __call__(self, phase, T, P):
+    def __call__(self, phase, T, P=None):
         if self.force_gas_critical_phase and T > self.Tc: phase = 'g'
         return getattr(self, phase)(T, P)
 

thermosteam/equilibrium/vle.py

@@ -1120,17 +1120,17 @@ def f(V):
                 self._T = thermal_condition.T = T
 
                 y0 = f(0.)
-                if y0 > self.H_hat_tol:
+                if y0 > 0:
                     self._T = thermal_condition.T = self.mixture.xsolve_T_at_HP(
                         self._phase_data, H, T, P
                     )
-                elif y0 < -self.H_hat_tol:
+                elif y0 < 0:
                     y1 = f(1.)
-                    if y1 < -self.H_hat_tol:
+                    if y1 < 0.:
                         self._T = thermal_condition.T = self.mixture.xsolve_T_at_HP(
                             self._phase_data, H, T, P
                         )
-                    elif y1 > self.H_hat_tol:
+                    elif y1 > 0.:
                         flx.IQ_interpolation(f,
                             0., 1., y0, y1, self._V, self.V_tol, self.H_hat_tol,
                             maxiter=self.maxiter, 

thermosteam/free_energy.py

@@ -22,96 +22,96 @@ def get_excess_energy(eos, T, P, free_energy, phase):
         return getattr(eos, name)
 
 @functor(var='H')
-def Enthalpy(T, Cn, T_ref, H_ref):
+def Enthalpy(T, P, Cn, T_ref, H_ref):
     return H_ref + Cn.T_dependent_property_integral(T_ref, T)
 
 @functor(var='S')
-def Entropy(T, Cn, T_ref, S0):
+def Entropy(T, P, Cn, T_ref, S0):
     return S0 + Cn.T_dependent_property_integral_over_T(T_ref, T)
 
 @functor(var='S')
 def EntropyGas(T, P, Cn, T_ref, P_ref, S0):
     return S0 + Cn.T_dependent_property_integral_over_T(T_ref, T) - R*log(P/P_ref)
 
 @functor(var='H.l')
-def Liquid_Enthalpy_Ref_Liquid(T, Cn_l, T_ref, H_ref):
+def Liquid_Enthalpy_Ref_Liquid(T, P, Cn_l, T_ref, H_ref):
     """Enthapy (kJ/kmol) disregarding pressure and assuming the specified phase."""
     return H_ref + Cn_l.T_dependent_property_integral(T_ref, T)
 
 @functor(var='H.l')
-def Liquid_Enthalpy_Ref_Gas(T, Cn_l, H_int_Tb_to_T_ref_g, Hvap_Tb, Tb, H_ref):
+def Liquid_Enthalpy_Ref_Gas(T, P, Cn_l, H_int_Tb_to_T_ref_g, Hvap_Tb, Tb, H_ref):
     return H_ref - H_int_Tb_to_T_ref_g - Hvap_Tb + Cn_l.T_dependent_property_integral(Tb, T)
 
 @functor(var='H.l')
 def Liquid_Enthalpy_Ref_Solid(T, Cn_l, H_int_T_ref_to_Tm_s, Hfus, Tm, H_ref):
     return H_ref + H_int_T_ref_to_Tm_s + Hfus + Cn_l.T_dependent_property_integral(Tm, T)
 
 @functor(var='H.s')
-def Solid_Enthalpy_Ref_Solid(T, Cn_s, T_ref, H_ref):
+def Solid_Enthalpy_Ref_Solid(T, P, Cn_s, T_ref, H_ref):
     return H_ref + Cn_s.T_dependent_property_integral(T_ref, T)
 
 @functor(var='H.s')
-def Solid_Enthalpy_Ref_Liquid(T, Cn_s, H_int_Tm_to_T_ref_l, Hfus, Tm, H_ref):
+def Solid_Enthalpy_Ref_Liquid(T, P, Cn_s, H_int_Tm_to_T_ref_l, Hfus, Tm, H_ref):
     return H_ref - H_int_Tm_to_T_ref_l - Hfus + Cn_s.T_dependent_property_integral(Tm, T)
 
 @functor(var='H.s')
-def Solid_Enthalpy_Ref_Gas(T, Cn_s, H_int_Tb_to_T_ref_g, Hvap_Tb,
+def Solid_Enthalpy_Ref_Gas(T, P, Cn_s, H_int_Tb_to_T_ref_g, Hvap_Tb,
                            H_int_Tm_to_Tb_l, Hfus, Tm, H_ref):
     return H_ref - H_int_Tb_to_T_ref_g - Hvap_Tb - H_int_Tm_to_Tb_l - Hfus + Cn_s.T_dependent_property_integral(Tm, T)
 
 @functor(var='H.g')
-def Gas_Enthalpy_Ref_Gas(T, Cn_g, T_ref, H_ref):
+def Gas_Enthalpy_Ref_Gas(T, P, Cn_g, T_ref, H_ref):
     return H_ref + Cn_g.T_dependent_property_integral(T_ref, T)
 
 @functor(var='H.g')
-def Gas_Enthalpy_Ref_Liquid(T, Cn_g, H_int_T_ref_to_Tb_l, Hvap_Tb, 
+def Gas_Enthalpy_Ref_Liquid(T, P, Cn_g, H_int_T_ref_to_Tb_l, Hvap_Tb, 
                             Tb, H_ref):
     return H_ref + H_int_T_ref_to_Tb_l + Hvap_Tb + Cn_g.T_dependent_property_integral(Tb, T)
 
 @functor(var='H.g')
-def Gas_Enthalpy_Ref_Solid(T, Cn_g, H_int_T_ref_to_Tm_s, Hfus,
+def Gas_Enthalpy_Ref_Solid(T, P, Cn_g, H_int_T_ref_to_Tm_s, Hfus,
                            H_int_Tm_to_Tb_l, Hvap_Tb, Tb, H_ref):
     return H_ref + H_int_T_ref_to_Tm_s + Hfus + H_int_Tm_to_Tb_l + Hvap_Tb + Cn_g.T_dependent_property_integral(Tb, T)
 
 
-EnthalpyRefLiquid = PhaseTFunctorBuilder('H',
+EnthalpyRefLiquid = PhaseTPFunctorBuilder('H',
                                         Solid_Enthalpy_Ref_Liquid.functor,
                                         Liquid_Enthalpy_Ref_Liquid.functor,
                                         Gas_Enthalpy_Ref_Liquid.functor)
 
-EnthalpyRefSolid = PhaseTFunctorBuilder('H',
+EnthalpyRefSolid = PhaseTPFunctorBuilder('H',
                                        Solid_Enthalpy_Ref_Solid.functor,
                                        Liquid_Enthalpy_Ref_Solid.functor,
                                        Gas_Enthalpy_Ref_Solid.functor)
 
-EnthalpyRefGas = PhaseTFunctorBuilder('H',
+EnthalpyRefGas = PhaseTPFunctorBuilder('H',
                                      Solid_Enthalpy_Ref_Gas.functor,
                                      Liquid_Enthalpy_Ref_Gas.functor,
                                      Gas_Enthalpy_Ref_Gas.functor)
 
 @functor(var='S.l')
-def Liquid_Entropy_Ref_Liquid(T, Cn_l, T_ref, S0):
+def Liquid_Entropy_Ref_Liquid(T, P, Cn_l, T_ref, S0):
     """Enthapy (kJ/kmol) disregarding pressure and assuming the specified phase."""
     return S0 + Cn_l.T_dependent_property_integral_over_T(T_ref, T)
 
 @functor(var='S.l')
-def Liquid_Entropy_Ref_Gas(T, Cn_l, S_int_Tb_to_T_ref_g, Svap_Tb, Tb, S0):
+def Liquid_Entropy_Ref_Gas(T, P, Cn_l, S_int_Tb_to_T_ref_g, Svap_Tb, Tb, S0):
     return S0 - S_int_Tb_to_T_ref_g - Svap_Tb + Cn_l.T_dependent_property_integral_over_T(Tb, T)
 
 @functor(var='S.l')
-def Liquid_Entropy_Ref_Solid(T, Cn_l, S_int_T_ref_to_Tm_s, Sfus, Tm, S0):
+def Liquid_Entropy_Ref_Solid(T, P, Cn_l, S_int_T_ref_to_Tm_s, Sfus, Tm, S0):
     return S0 + S_int_T_ref_to_Tm_s + Sfus + Cn_l.T_dependent_property_integral_over_T(Tm, T)
 
 @functor(var='S.s')
-def Solid_Entropy_Ref_Solid(T, Cn_s, T_ref, S0):
+def Solid_Entropy_Ref_Solid(T, P, Cn_s, T_ref, S0):
     return S0 + Cn_s.T_dependent_property_integral_over_T(T_ref, T)
 
 @functor(var='S.s')
-def Solid_Entropy_Ref_Liquid(T, Cn_s, S_int_Tm_to_T_ref_l, Sfus, Tm, S0):
+def Solid_Entropy_Ref_Liquid(T, P, Cn_s, S_int_Tm_to_T_ref_l, Sfus, Tm, S0):
     return S0 - S_int_Tm_to_T_ref_l - Sfus + Cn_s.T_dependent_property_integral_over_T(Tm, T)
 
 @functor(var='S.s')
-def Solid_Entropy_Ref_Gas(T, Cn_s, S_int_Tb_to_T_ref_g, Svap_Tb, 
+def Solid_Entropy_Ref_Gas(T, P, Cn_s, S_int_Tb_to_T_ref_g, Svap_Tb, 
                           S_int_Tm_to_Tb_l, Sfus, Tm, S0):
     return S0 - S_int_Tb_to_T_ref_g - Svap_Tb - S_int_Tm_to_Tb_l - Sfus + Cn_s.T_dependent_property_integral_over_T(Tm, T)
 
@@ -245,7 +245,7 @@ def Excess_Gas_Entropy_Ref_Liquid(T, P, eos, S_dep_T_ref_Pb,
     return S_dep_T_ref_Pb - S_dep_ref_l + S_dep_g - S_dep_Tb_Pb_g
 
 @functor(var='S.g')
-def Excess_Gas_Entropy_Ref_Solid(T):
+def Excess_Gas_Entropy_Ref_Solid(T, P):
     return 0
 
 ExcessEntropyRefLiquid = PhaseTPFunctorBuilder('S_excess',

thermosteam/mixture/mixture.py

@@ -34,7 +34,7 @@ def create_mixture_model(chemicals, var, Model):
         obj = getfield(chemical, var)
         if isa(obj, PhaseHandle):
             phase_handle = obj
-        elif var in ('Cn', 'H'):
+        elif var == 'Cn':
             phase_handle = MockPhaseTHandle(var, obj)
         else:
             phase_handle = MockPhaseTPHandle(var, obj)
@@ -226,7 +226,7 @@ def from_chemicals(cls, chemicals,
                 MWs = chemical_data_array(chemicals, 'MW')
             getfield = getattr
             Cn =  create_mixture_model(chemicals, 'Cn', IdealTMixtureModel)
-            H =  create_mixture_model(chemicals, 'H', IdealTMixtureModel)
+            H =  create_mixture_model(chemicals, 'H', IdealTPMixtureModel)
             S = create_mixture_model(chemicals, 'S', IdealEntropyModel)
             H_excess = create_mixture_model(chemicals, 'H_excess', IdealTPMixtureModel)
             S_excess = create_mixture_model(chemicals, 'S_excess', IdealTPMixtureModel)
@@ -321,7 +321,7 @@ def get_property(self, name, units, *args, **kwargs):
 
     def H(self, phase, mol, T, P):
         """Return enthalpy [J/mol]."""
-        H = self._H(phase, mol, T)
+        H = self._H(phase, mol, T, P)
         if self.include_excess_energies: H += self._H_excess(phase, mol, T, P)
         return H
 
@@ -396,7 +396,7 @@ def xH(self, phase_mol, T, P):
         """Multi-phase mixture enthalpy [J/mol]."""
         H = self._H
         phase_mol = tuple(phase_mol)
-        H_total = sum([H(phase, mol, T) for phase, mol in phase_mol])
+        H_total = sum([H(phase, mol, T, P) for phase, mol in phase_mol])
         if self.include_excess_energies:
             H_excess = self._H_excess
             H_total += sum([H_excess(phase, mol, T, P) for phase, mol in phase_mol])