implement ideal-gas poyinting correction; minor name chage

thermosteam/_chemical.py

@@ -419,7 +419,7 @@ class Chemical:
     Phase dependent properties have attributes with model handles for each phase:
 
     >>> Water.V
-    <PhaseTPHandle(phase, T, P) -> V [m^3/mol]>
+    <PhaseTPHandle(phase, T, P=None) -> V [m^3/mol]>
     >>> # Water.V.l
     >>> # VolumeLiquid(CASRN="7732-18-5", MW=18.01528, Tb=373.124, Tc=647.14, Pc=22048320.0, Vc=5.6-05, Zc=0.2294727, omega=0.344, dipole=1.85, Psat=VaporPressure(CASRN="7732-18-5", Tb=373.124, Tc=647.14, Pc=22048320.0, omega=0.344, extrapolation="AntoineAB|DIPPR101_ABC", method="WAGNER_MCGARRY"), eos=[PR(Tc=647.14, Pc=22048320.0, omega=0.344, T=298.15, P=101325.0)], extrapolation="constant", method="VDI_PPDS", method_P="COSTALD_COMPRESSED", tabular_extrapolation_permitted=True)
     
@@ -442,7 +442,7 @@ class Chemical:
 
     Choose what model to use through the `method` attribute:
     
-    >>> list(sorted(Water.Cn.l.all_methods))
+    >>> list(sorted(Water.Cn.l.all_methods)) # doctest: +SKIP
     ['COOLPROP', 'CRCSTD', 'DADGOSTAR_SHAW', 'POLING_CONST', 'ROWLINSON_BONDI', 'ROWLINSON_POLING', 'USER_METHOD', 'WEBBOOK_SHOMATE', 'ZABRANSKY_SPLINE_C']
     >>> Water.Cn.l.method = 'ZABRANSKY_SPLINE_C'
     
@@ -2072,7 +2072,7 @@ def at_state(self, phase, copy=False):
         >>> from thermosteam import Chemical
         >>> N2 = Chemical('N2')
         >>> N2.at_state(phase='g')
-        >>> N2.H(298.15) # No longer a function of phase
+        >>> N2.H(298.15, None) # No longer a function of phase
         0.0
         """
         if copy:

thermosteam/_settings.py

@@ -64,7 +64,7 @@ class ProcessSettings:
         ),
         Gamma=DortmundActivityCoefficients,
         Phi=IdealFugacityCoefficients,
-        PCF=IdealPoyintingCorrectionFactors
+        PCF=MockPoyintingCorrectionFactors
     )
 
     Access defined chemicals:

thermosteam/_thermo.py

@@ -55,7 +55,7 @@ class Thermo:
         ),
         Gamma=DortmundActivityCoefficients,
         Phi=IdealFugacityCoefficients,
-        PCF=IdealPoyintingCorrectionFactors
+        PCF=MockPoyintingCorrectionFactors
     )
     
     Note that the Dortmund-UNIFAC is the default activity coefficient model. 
@@ -108,7 +108,7 @@ class Thermo:
         ),
         Gamma=DortmundActivityCoefficients,
         Phi=IdealFugacityCoefficients,
-        PCF=IdealPoyintingCorrectionFactors
+        PCF=MockPoyintingCorrectionFactors
     )
     
     Attributes
@@ -130,7 +130,7 @@ class Thermo:
     def __init__(self, chemicals, mixture=None,
                  Gamma=eq.DortmundActivityCoefficients,
                  Phi=eq.IdealFugacityCoefficients,
-                 PCF=eq.IdealPoyintingCorrectionFactors,
+                 PCF=eq.MockPoyintingCorrectionFactors,
                  cache=None,
                  skip_checks=False):
         if not isinstance(chemicals, Chemicals): chemicals = Chemicals(chemicals, cache)
@@ -274,7 +274,7 @@ class IdealThermo:
 
     Gamma = eq.IdealActivityCoefficients
     Phi = eq.IdealFugacityCoefficients
-    PCF = eq.IdealPoyintingCorrectionFactors
+    PCF = eq.MockPoyintingCorrectionFactors
     as_chemical = Thermo.as_chemical
     __enter__ = Thermo.__enter__
     __exit__ = Thermo.__exit__

thermosteam/base/phase_handle.py

@@ -127,7 +127,7 @@ def __call__(self, phase, T, P=None):
 class MockPhaseTPHandle(MockPhaseHandle):
     __slots__ = ()
 
-    def __call__(self, phase, T, P):
+    def __call__(self, phase, T, P=None):
         return self.model(T, P)
 
 # %% Builders

thermosteam/equilibrium/bubble_point.py

@@ -114,21 +114,22 @@ def __init__(self, chemicals=(), thermo=None):
 
     def _T_error(self, T, P, z_over_P, z_norm, y):
         if T <= 0: raise InfeasibleRegion('negative temperature')
+        Psats = np.array([i(T) for i in self.Psats], dtype=float)
         y_phi =  (z_over_P
-                  * np.array([i(T) for i in self.Psats])
+                  * Psats
                   * self.gamma(z_norm, T) 
-                  * self.pcf(T, P))
+                  * self.pcf(T, P, Psats))
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
 
-    def _P_error(self, P, T, z_Psat_gamma, y):
+    def _P_error(self, P, T, z_Psat_gamma, Psats, y):
         if P <= 0: raise InfeasibleRegion('negative pressure')
-        y_phi = z_Psat_gamma * self.pcf(T, P) / P
+        y_phi = z_Psat_gamma * self.pcf(T, P, Psats) / P
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
 
     def _T_error_ideal(self, T, z_over_P, y):
-        y[:] = z_over_P * np.array([i(T) for i in self.Psats])
+        y[:] = z_over_P * np.array([i(T) for i in self.Psats], dtype=float)
         return 1 - y.sum()
 
     def _Ty_ideal(self, z_over_P):
@@ -263,7 +264,7 @@ def solve_Py(self, z, T):
             z_Psat_gamma = z * Psat * self.gamma(z_norm, T)
             f = self._P_error
             P_guess, y = self._Py_ideal(z_Psat_gamma)
-            args = (T, z_Psat_gamma, y)
+            args = (T, z_Psat_gamma, Psat, y)
             try:
                 P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9,
                                       args, checkiter=False)

thermosteam/equilibrium/dew_point.py

@@ -21,30 +21,29 @@
 # %% Solvers
 
 # @njit(cache=True)
-def x_iter(x, x_gamma_poyinting, T, P, f_gamma, gamma_args, f_pcf, pcf_args):
+def x_iter(x, x_gamma, T, P, f_gamma, gamma_args):
     # Add back trace amounts for activity coefficients at infinite dilution
     mask = x < 1e-32
     x[mask] = 1e-32
     x = fn.normalize(x)
     gamma = f_gamma(x, T, *gamma_args)
-    pcf = f_pcf(T, P, *pcf_args)
-    denominator = gamma * pcf
+    denominator = gamma
     try:
-        x = x_gamma_poyinting / denominator
+        x = x_gamma / denominator
     except: 
         raise Exception('liquid phase composition is infeasible')
     if (np.abs(x) > 1e16).any():
         raise Exception('liquid phase composition is infeasible')
     return x
 
 # @njit(cache=True)
-def solve_x(x_guess, x_gamma_poyinting, T, P, f_gamma, gamma_args, f_pcf, pcf_args):
-    args = (x_gamma_poyinting, T, P, f_gamma, gamma_args, f_pcf, pcf_args)
+def solve_x(x_guess, x_gamma, T, P, f_gamma, gamma_args):
+    args = (x_gamma, T, P, f_gamma, gamma_args)
     try:
         x = flx.wegstein(x_iter, x_guess, 1e-10, args=args, checkiter=False,
                          checkconvergence=False, convergenceiter=3)
     except:
-        return x_gamma_poyinting
+        return x_gamma
     return x
 
 # %% Dew point values container
@@ -125,18 +124,18 @@ def __init__(self, chemicals=(), thermo=None):
             self.chemicals = chemicals
             cached[key] = self
 
-    def _solve_x(self, x_gamma_pcf, T, P, x):
+    def _solve_x(self, x_gamma, T, P, x):
         gamma = self.gamma
-        pcf = self.pcf
-        return solve_x(x, x_gamma_pcf, T, P, gamma.f, gamma.args, pcf.f, pcf.args)
+        return solve_x(x, x_gamma, T, P, gamma.f, gamma.args)
 
     def _T_error(self, T, P, z_norm, zP, x):
         if T <= 0: raise InfeasibleRegion('negative temperature')
         Psats = np.array([i(T) for i in self.Psats])
         Psats[Psats < 1e-16] = 1e-16 # Prevent floating point error
         phi = self.phi(z_norm, T, P)
-        x_gamma_pcf = phi * zP / Psats
-        x[:] = self._solve_x(x_gamma_pcf, T, P, x)
+        pcf = self.pcf(T, P, Psats)
+        x_gamma = phi * zP / Psats / pcf
+        x[:] = self._solve_x(x_gamma, T, P, x)
         return 1 - x.sum()
 
     def _T_error_ideal(self, T, zP, x):
@@ -145,10 +144,10 @@ def _T_error_ideal(self, T, zP, x):
         x[:] = zP / Psats
         return 1 - x.sum()
 
-    def _P_error(self, P, T, z_norm, z_over_Psats, x):
+    def _P_error(self, P, T, z_norm, z_over_Psats, Psats, x):
         if P <= 0: raise InfeasibleRegion('negative pressure')
-        x_gamma_pcf = z_over_Psats * P * self.phi(z_norm, T, P)
-        x[:] = self._solve_x(x_gamma_pcf, T, P, x)
+        x_gamma = z_over_Psats * P * self.phi(z_norm, T, P) / self.pcf(T, P, Psats)
+        x[:] = self._solve_x(x_gamma, T, P, x)
         return 1 - x.sum()
 
     def _Tx_ideal(self, zP):
@@ -280,7 +279,7 @@ def solve_Px(self, z, T):
             Psats = np.array([i(T) for i in self.Psats], dtype=float)
             z_over_Psats = z/Psats
             P_guess, x = self._Px_ideal(z_over_Psats)
-            args = (T, z_norm, z_over_Psats, x)
+            args = (T, z_norm, z_over_Psats, Psats, x)
             f = self._P_error
             try:
                 P = flx.aitken_secant(f, P_guess, P_guess-10, self.P_tol, 5e-12, args,

thermosteam/equilibrium/fugacities.py

@@ -51,7 +51,8 @@ def __init__(self, chemicals, thermo=None):
         self.pcf = thermo.PCF(chemicals)
 
     def __call__(self, x, T, P=101325.):
-        return x * self.gamma(x, T) * self.pcf(T, P) * np.array([i.Psat(T) for i in self.chemicals], dtype=float)
+        Psats = np.array([i.Psat(T) for i in self.chemicals], dtype=float)
+        return x * self.gamma(x, T) * self.pcf(T, P, Psats) * Psats
 
     def __repr__(self):
         chemicals = ", ".join([i.ID for i in self.chemicals])

thermosteam/equilibrium/ideal.py

@@ -6,7 +6,7 @@
 """
 from numba import njit
 
-__all__ = ('ideal_coefficient',)
+__all__ = ('ideal',)
 
 def ideal(cls):
     cls.f = ideal_coefficient

thermosteam/equilibrium/poyinting_correction_factors.py

@@ -8,47 +8,79 @@
 """
 """
 __all__ = ('PoyintingCorrectionFactors',
-           'IdealPoyintingCorrectionFactors')
-
-from .ideal import ideal
+           'MockPoyintingCorrectionFactors',
+           'IdealGasPoyintingCorrectionFactors')
+from ..constants import R
+from numba import njit
+import numpy as np
 
 class PoyintingCorrectionFactors:
-    """Abstract class for the estimation of Poyinting correction factors. Non-abstract subclasses should implement the following methods:
+    """Abstract class for the estimation of Poyinting correction factors.
+    Non-abstract subclasses should implement the following methods:
         
     __init__(self, chemicals: Iterable[:class:`~thermosteam.Chemical`]):
-        Should use pure component data from chemicals to setup future calculations of Poyinting correction factors.
+        Should use pure component data from chemicals to setup future 
+        calculations of Poyinting correction factors.
     
     __call__(self, T: float, P: float):
-        Should accept the temperature `T` (in Kelvin) and pressure `P` (in Pascal), and return an array of Poyinting correction factors. 
+        Should accept the temperature `T` (in Kelvin) and pressure `P` (in Pascal),
+        and return an array of Poyinting correction factors. 
         
     """
     __slots__ = ()
 
+    @property
+    def chemicals(self):
+        return self._chemicals
+    @chemicals.setter
+    def chemicals(self, chemicals):
+        self._chemicals = tuple(chemicals)
+
     def __init__(self, chemicals):
         self.chemicals = chemicals
 
+    def __reduce__(self):
+        return type(self), (self.chemicals,)
+
     def __repr__(self):
         chemicals = ", ".join([i.ID for i in self.chemicals])
         return f"<{type(self).__name__}([{chemicals}])>"
 
-@ideal
-class IdealPoyintingCorrectionFactors(PoyintingCorrectionFactors):
-    """Create an IdealPoyintingCorrectionFactor object that estimates all poyinting correction factors to be 1 when called with composition and temperature (K).
+
+class MockPoyintingCorrectionFactors(PoyintingCorrectionFactors):
+    """Create a MockPoyintingCorrectionFactors object 
+    that estimates all poyinting correction factors to be 1 when
+    called with a temperature (K) and pressure (Pa).
     
     Parameters
     ----------
     
     chemicals : Iterable[:class:`~thermosteam.Chemical`]
     
     """
-    __slots__ = ('_chemicals')
+    __slots__ = ('_chemicals',)
 
-    @property
-    def chemicals(self):
-        return self._chemicals
-    @chemicals.setter
-    def chemicals(self, chemicals):
-        self._chemicals = tuple(chemicals)
-
-    def __call__(self, T, P):
+    def __call__(self, T, P, Psats):
         return 1.
+
+
+@njit(cache=True)
+def ideal_gas_poyinting_correction_factors(T, P, vls, Psats):
+    return np.exp(vls / (R * T) * (P - Psats))
+
+class IdealGasPoyintingCorrectionFactors(PoyintingCorrectionFactors):
+    """Create an IdealGasPoyintingCorrectionFactors object that estimates 
+    poyinting correction factors assuming ideal gas when called with 
+    a temperature (K) and pressure (Pa).
+    
+    Parameters
+    ----------
+    
+    chemicals : Iterable[:class:`~thermosteam.Chemical`]
+    
+    """
+    __slots__ = ('_chemicals',)
+
+    def __call__(self, T, P, Psats):
+        vls = np.array([i.V('l', T, P) for i in self._chemicals])
+        return ideal_gas_poyinting_correction_factors(T, P, vls, Psats)