improvements to LLE; workaround for bug (needs future fix)

Author: yoelcortes

thermosteam/equilibrium/bubble_point.py

@@ -39,22 +39,51 @@ def solve_y(y_phi, phi, T, P, y_guess):
 # %% Bubble point values container
 
 class BubblePointValues:
-    __slots__ = ('T', 'P', 'IDs', 'z', 'y')
+    __slots__ = ('T', 'P', 'IDs', 'z', 'y', 'α')
 
-    def __init__(self, T, P, IDs, z, y):
+    def __init__(self, T, P, IDs, z, y, α=None):
         self.T = T
         self.P = P
         self.IDs = IDs
         self.z = z
         self.y = y
-        
+        self.α = α
+    
     @property
     def x(self): return self.z
 
     @property
     def K(self):
         return self.y / self.x
 
+    @property
+    def xα(self):
+        return self.α / self.α.sum()
+    
+    @property
+    def xβ(self):
+        return self.β / self.β.sum()
+    
+    @property
+    def β(self):
+        return self.z - self.α
+    
+    @property
+    def ϕ(self):
+        return self.α.sum()
+    
+    @property
+    def KL(self):
+        return self.xα / self.xβ
+    
+    @property
+    def Kα(self):
+        return self.y / self.xα
+    
+    @property
+    def Kβ(self):
+        return self.y / self.xβ
+    
     def __repr__(self):
         return f"{type(self).__name__}(T={self.T:.2f}, P={self.P:.0f}, IDs={self.IDs}, z={self.z}, y={self.y})"
 
@@ -157,11 +186,14 @@ def _T_error(self, T, P, z_over_P, z_norm, y):
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
 
-    def _T_error_lle(self, T, P, z, y, x):
+    def _T_error_lle(self, T, P, z, y, x, α):
         if T <= 0: raise InfeasibleRegion('negative temperature')
         Psats = np.array([i(T) for i in self.Psats], dtype=float)
-        mol = solve_lle_mol(self.gamma, z, T, P, sample=x)
-        x[:] = mol / mol.sum() if mol.any() else z
+        α[:] = solve_lle_mol(self.gamma, z, T, P, sample=x)
+        if α.any(): 
+            x[:] = α / α.sum()
+        else:
+            x = z
         y_phi =  (x / P
                   * Psats
                   * self.gamma(x, T, P) 
@@ -181,10 +213,13 @@ def _P_error_dep(self, P, T, z, Psats, z_Psats, y):
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
 
-    def _P_error_lle(self, P, T, z, Psats, y, x):
+    def _P_error_lle(self, P, T, z, Psats, y, x, α):
         if P <= 0: raise InfeasibleRegion('negative pressure')
-        mol = solve_lle_mol(self.gamma, z, T, P, sample=x)
-        x[:] = mol / mol.sum() if mol.any() else z
+        α[:] = solve_lle_mol(self.gamma, z, T, P, sample=x)
+        if α.any(): 
+            x[:] = α / α.sum()
+        else:
+            x = z
         y_phi = Psats * x * self.gamma(x, T, P) * self.pcf(T, P, Psats) / P
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
@@ -242,17 +277,17 @@ def __call__(self, z, *, T=None, P=None, liquid_conversion=None, lle=False):
         z = np.asarray(z, float)
         if T:
             if P: raise ValueError("may specify either T or P, not both")
-            P, *args = self.solve_Py(z, T, liquid_conversion, lle)
+            P, *args = self.solve_Py(z, T, liquid_conversion, lle, full=True)
         elif P:
-            T, *args = self.solve_Ty(z, P, liquid_conversion, lle)
+            T, *args = self.solve_Ty(z, P, liquid_conversion, lle, full=True)
         else:
             raise ValueError("must specify either T or P")
         if liquid_conversion:
             return ReactiveBubblePointValues(T, P, self.IDs, z, *args)
         else:
             return BubblePointValues(T, P, self.IDs, z, *args)
 
-    def solve_Ty(self, z, P, liquid_conversion=None, lle=False):
+    def solve_Ty(self, z, P, liquid_conversion=None, lle=False, full=False):
         """
         Bubble point at given composition and pressure.
 
@@ -295,10 +330,13 @@ def solve_Ty(self, z, P, liquid_conversion=None, lle=False):
             z_over_P = z / P
             T_guess, y = self._Ty_ideal(z_over_P)
             if lle:
+                x = y.copy()
+                α = x.copy()
                 f = self._T_error_lle
-                args = (P, z, y, y.copy())
+                args = (P, z, y, x, α)
             else:
                 f = self._T_error
+                α= None
                 args = (P, z_over_P, z, y)
             try:
                 T = flx.aitken_secant(f, T_guess, T_guess + 1e-3,
@@ -313,7 +351,8 @@ def solve_Ty(self, z, P, liquid_conversion=None, lle=False):
                                          T_guess, self.T_tol, 5e-12, args, 
                                          checkiter=False, checkbounds=False, 
                                          maxiter=self.maxiter)
-            return T, fn.normalize(y)
+            y = fn.normalize(y)
+            return (T, y, α) if full else (T, y)
         else:
             f = self._T_error_reactive
             z_norm = z / z.sum()
@@ -336,7 +375,7 @@ def solve_Ty(self, z, P, liquid_conversion=None, lle=False):
                                          maxiter=self.maxiter)
             return T, dz, fn.normalize(y), x
 
-    def solve_Py(self, z, T, liquid_conversion=None, lle=False):
+    def solve_Py(self, z, T, liquid_conversion=None, lle=False, full=False):
         """
         Bubble point at given composition and temperature.
 
@@ -379,6 +418,7 @@ def solve_Py(self, z, T, liquid_conversion=None, lle=False):
             elif T < self.Tmin: T = self.Tmin
             Psats = np.array([i(T) for i in self.Psats])
             z = z / z.sum()
+            α = None
             if self.gamma.P_dependent:
                 f = self._P_error_dep
                 z_Psats = z * Psats
@@ -389,7 +429,9 @@ def solve_Py(self, z, T, liquid_conversion=None, lle=False):
                     z_Psat_gamma = z * Psats * self.gamma(z, T, 101325)
                     P_guess, y = self._Py_ideal(z_Psat_gamma)
                     f = self._P_error_lle
-                    args = (P, T, z, Psats, y, y.copy())
+                    x = y.copy()
+                    α = x.copy()
+                    args = (P, T, z, Psats, y, x, α)
                 else:
                     z_Psat_gamma = z * Psats * self.gamma(z, T, 101325)
                     P_guess, y = self._Py_ideal(z_Psat_gamma)
@@ -406,7 +448,8 @@ def solve_Py(self, z, T, liquid_conversion=None, lle=False):
                                          P_guess, self.P_tol, 5e-12, args,
                                          checkiter=False, checkbounds=False, 
                                          maxiter=self.maxiter)
-            return P, fn.normalize(y)
+            y = fn.normalize(y)
+            return (P, y, α) if full else P, y
         else:
             f = self._P_error_reactive
             z_norm = z / z.sum()

thermosteam/equilibrium/lle.py

@@ -135,11 +135,14 @@ def solve_lle_mol(gamma, z, T, P, sample=None, phi=1):
     else:
         return z
     K[K <= 0] = 1e-9
-    mol = pseudo_equilibrium(
-        K, phi, z, T, n, gamma.f, gamma.args, 
-        LLE.pseudo_equilibrium_inner_loop_options,
-        LLE.pseudo_equilibrium_outer_loop_options,
-    )
+    try:
+        mol = pseudo_equilibrium(
+            K, phi, z, T, n, gamma.f, gamma.args, 
+            LLE.pseudo_equilibrium_inner_loop_options,
+            LLE.pseudo_equilibrium_outer_loop_options,
+        )
+    except:
+        breakpoint()
     if mol.any():
         return mol / mol.sum()
     else:

thermosteam/equilibrium/plot_equilibrium.py

@@ -60,7 +60,7 @@ def as_thermo(thermo, chemicals): # pragma: no cover
 # %% Plot functions
 
 def plot_vle_binary_phase_envelope(chemicals, T=None, P=None, vc=None, lc=None, thermo=None,
-                                   yticks=None, line_styles=None, N=None): # pragma: no cover
+                                   yticks=None, line_styles=None, N=None, lle=False): # pragma: no cover
     """
     Plot the binary phase envelope of two chemicals at a given temperature or pressure.
 
@@ -100,15 +100,15 @@ def plot_vle_binary_phase_envelope(chemicals, T=None, P=None, vc=None, lc=None,
         Xs = P
         X_units = 'Pa'
         assert not T, "must pass either T or P, but not both"
-        bpss = [[BP(z, P=i) for z in zs] for i in P]
+        bpss = [[BP(z, P=i, lle=lle) for z in zs] for i in P]
         mss = [[bp.T for bp in bps] for bps in bpss]
         ylabel = 'Temperature [K]'
     elif T:
         if not isinstance(T, Iterable): T = [T]
         assert not P, "must pass either T or P, but not both"
         Xs = T
         X_units = 'K'
-        bpss = [[BP(z, T=i) for z in zs] for i in T]
+        bpss = [[BP(z, T=i, lle=lle) for z in zs] for i in T]
         mss = [[bp.P for bp in bps] for bps in bpss]
         ylabel = 'Pressure [Pa]'
     else:
@@ -118,7 +118,7 @@ def plot_vle_binary_phase_envelope(chemicals, T=None, P=None, vc=None, lc=None,
     for X, ms, bps, ls in zip(Xs, mss, bpss, line_styles):
         ys_a = np.array([bp.y[0] for bp in bps])
         ms_liq = ms.copy()
-        ms_gas = interpolate.interp1d(ys_a, ms, bounds_error=False, kind='slinear')(zs_a)
+        # ms_gas = interpolate.interp1d(ys_a, ms, bounds_error=False, kind='slinear')(zs_a)
         # if P:
         #     azeotrope = ms_liq > ms_gas
         # elif T:
@@ -134,7 +134,7 @@ def plot_vle_binary_phase_envelope(chemicals, T=None, P=None, vc=None, lc=None,
         #     azeotrope[mid] = False
         #     ms_gas[mid] = ms_liq[mid]
         #     ms_gas = interpolate.interp1d(zs_a[~azeotrope], ms_gas[~azeotrope], bounds_error=False, kind='slinear')(zs_a)
-        plt.plot(zs_a, ms_gas, c=vc if vc is not None else colors.red.RGBn, label=f'vapor [{int(X)} {X_units}]', ls=ls)
+        plt.plot(ys_a, ms, c=vc if vc is not None else colors.red.RGBn, label=f'vapor [{int(X)} {X_units}]', ls=ls)
         plt.plot(zs_a, ms_liq, c=lc if lc is not None else colors.blue.RGBn, label=f'liquid [{int(X)} {X_units}]', ls=ls)
     plt.xlim([0, 1])
     plt.ylim([mss.min(), mss.max()])

thermosteam/equilibrium/tangential_plane_stability.py

@@ -22,7 +22,7 @@ class StabilityReport(NamedTuple):
     sample_unstable: bool = False
 
 def edge_points_simplex_masked(z: np.ndarray,
-                               points_per_edge: int = 8,
+                               points_per_edge: int = 12,
                                epsilon: float = 1e-3,
                                min_active: int = 2) -> np.ndarray:
     """
@@ -121,14 +121,6 @@ def objective(w, T, P, logfz, softmax=False):
             w /= w.sum()
             best_result = w
             best_val = value
-        unstable = best_val < 0
-        if unstable:
-            return StabilityReport(
-                unstable=unstable,
-                candidate=sample,
-                tpd=best_val,
-                sample_unstable=True,
-            )
     w = z * MW
     w /= w.sum()
     samples = edge_points_simplex_masked(w)
@@ -154,7 +146,7 @@ def objective(w, T, P, logfz, softmax=False):
                 best_result = w
                 best_val = value
     return StabilityReport(
-        unstable=(best_val < 0),
+        unstable=best_val < 0 and np.abs(best_result - z).sum() > 1e-9,
         candidate=best_result,
         tpd=best_val
     )
@@ -202,14 +194,6 @@ def __call__(self, z, T, P, reference_phase='l', potential_phase='L', sample=Non
                 w /= w.sum()
                 best_result = w
                 best_val = value
-            unstable = best_val < 0
-            if unstable:
-                return StabilityReport(
-                    unstable=unstable,
-                    candidate=sample,
-                    tpd=best_val,
-                    sample_unstable=True,
-                )
         MW = self.MW
         w = z * MW
         w /= w.sum()
@@ -236,7 +220,7 @@ def __call__(self, z, T, P, reference_phase='l', potential_phase='L', sample=Non
                     best_result = w
                     best_val = value
         return StabilityReport(
-            unstable=(best_val < 0),
+            unstable=best_val < 0 and np.abs(best_result - z).sum() > 1e-9,
             candidate=best_result,
             tpd=best_val
         )

thermosteam/indexer.py

@@ -1423,20 +1423,26 @@ def by_mass(self):
     """Return a ChemicalMassFlowIndexer that references this object's molar data."""
     try:
         mass = self._data_cache['mass']
-    except:
-        chemicals = self.chemicals
-        dct = MassFlowDict(self.data.dct, chemicals.MW)
-        data_cache = self._data_cache
-        data_cache['mass'] = mass = \
-        ChemicalMassFlowIndexer.from_data(
-            SparseVector.from_dict(
-                dct,
-                chemicals.size
-            ),
-            self, 
-            chemicals,
-            False,
-        )
+        # TODO: understand where is this error coming from.
+        # Create test case trying to break it.
+        if mass.data.dct.dct is self.data.dct:
+            return mass
+        else:
+            self._data_cache.clear()
+    except: pass
+    chemicals = self.chemicals
+    dct = MassFlowDict(self.data.dct, chemicals.MW)
+    data_cache = self._data_cache
+    data_cache['mass'] = mass = \
+    ChemicalMassFlowIndexer.from_data(
+        SparseVector.from_dict(
+            dct,
+            chemicals.size
+        ),
+        self, 
+        chemicals,
+        False,
+    )
     return mass
 ChemicalMolarFlowIndexer.by_mass = by_mass
 

thermosteam/units_of_measure.py

@@ -114,6 +114,9 @@ def format_units(units, ends='$', mathrm=True):
     all_numerators = []
     all_denominators = []
     unprocessed_numerators, *unprocessed_denominators = units.split("/")
+    unprocessed_numerators = unprocessed_numerators.strip(' ')
+    unprocessed_numerators = unprocessed_numerators.replace(' \cdot ', '*')
+    unprocessed_numerators = unprocessed_numerators.replace(' * ', '*')
     unprocessed_numerators = unprocessed_numerators.replace(' ', '*')
     unprocessed_numerators = unprocessed_numerators.replace('*%', '\ %')
     unprocessed_numerators = unprocessed_numerators.replace('%*', '%\ ')