minor improvements

thermosteam/_chemicals.py

@@ -489,6 +489,14 @@ def chemical_groups(self) -> frozenset[str]:
         """All defined chemical groups."""
         return frozenset(self._group_mol_compositions)
 
+    def chemical_group_members(self, name) -> list[str]:
+        index = self._index
+        if name in index:
+            IDs = self.IDs
+            return [IDs[i] for i in index[name]]
+        else:
+            raise ValueError(f"chemical group '{name}' does not exist")
+
     def refresh_constants(self):
         """
         Refresh constant arrays according to their chemical values,
@@ -536,7 +544,8 @@ def _compile(self, chemicals, skip_checks=False):
             raise RuntimeError(
                 f"{chemical} is missing key thermodynamic properties ({missing}); "
                 "use the `<Chemical>.get_missing_properties()` to check "
-                "all missing properties")
+                "all missing properties"
+            )
         IDs = tuple_([i.ID for i in chemicals])
         CAS = tuple_([i.CAS for i in chemicals])
         size = len(IDs)

thermosteam/_stream.py

@@ -449,7 +449,6 @@ def _update_energy_departure_coefficient(self, coefficients):
         coefficients[key] = value
 
     def _update_material_flows(self, value, index=None):
-        value[value < 0] = 0
         if index is None:
             self.mol[:] = value
         else:

thermosteam/equilibrium/lle.py

@@ -195,7 +195,7 @@ def __init__(self, imol=None, thermal_condition=None, thermo=None,
         self._K = None
         self._phi = None
 
-    def __call__(self, T, P=None, top_chemical=None, update=True, use_cache=True):
+    def __call__(self, T, P=None, top_chemical=None, update=True, use_cache=True, single_loop=False):
         """
         Perform liquid-liquid equilibrium.
 
@@ -242,7 +242,7 @@ def __call__(self, T, P=None, top_chemical=None, update=True, use_cache=True):
                 if self._lle_chemicals != lle_chemicals: 
                     self._K = None
                     self._phi = None
-                mol_L = self.solve_lle_liquid_mol(mol, T, lle_chemicals)
+                mol_L = self.solve_lle_liquid_mol(mol, T, lle_chemicals, single_loop)
                 mol_l = mol - mol_L
             if top_chemical:
                 MW = self.chemicals.MW[index]
@@ -302,7 +302,7 @@ def __call__(self, T, P=None, top_chemical=None, update=True, use_cache=True):
                 phi = 0.
             return lle_chemicals, K, phi
 
-    def solve_lle_liquid_mol(self, mol, T, lle_chemicals):
+    def solve_lle_liquid_mol(self, mol, T, lle_chemicals, single_loop):
         gamma = self.thermo.Gamma(lle_chemicals)
         indices = np.argsort(mol * np.array([i.MW for i in lle_chemicals]))
         method = self.method
@@ -324,11 +324,43 @@ def solve_lle_liquid_mol(self, mol, T, lle_chemicals):
                 y /= y.sum()
                 K = gamma(y, T) / gamma(x, T)
                 phi = 0.5
-            return pseudo_equilibrium(
-                K, phi, mol, T, n, gamma.f, gamma.args, 
-                self.pseudo_equilibrium_inner_loop_options,
-                self.pseudo_equilibrium_outer_loop_options,
-            )
+            if single_loop:
+                z = mol
+                f_gamma = gamma.f
+                gamma_args = gamma.args
+                phi = phase_fraction(z, K, phi)
+                try:
+                    x = z/(1. + phi * (K - 1.))
+                except:
+                    x = np.ones(n)
+                x /= x.sum()
+                y = K * x
+                Kgammayphi = np.zeros(2*n + 1)
+                Kgammayphi[:n] = K
+                Kgammayphi[n:-1] = f_gamma(y, T, *gamma_args)
+                Kgammayphi[-1] = phi
+                Kgammay = Kgammayphi[:-1]
+                phi = Kgammayphi[-1]
+                args=(z, T, n, f_gamma, gamma_args)
+                Kgammay = flx.fixed_point(
+                    psuedo_equilibrium_inner_loop, Kgammay, 
+                    args=(*args, phi), **self.pseudo_equilibrium_inner_loop_options,
+                )
+                K = Kgammay[:n]
+                try:
+                    phi = phase_fraction(z, K, phi)
+                except (ZeroDivisionError, FloatingPointError):
+                    return z
+                if np.isnan(phi): return z
+                if phi > 1: phi = 1 - 1e-16
+                if phi < 0: phi = 1e-16
+                return z/(1. + phi * (K - 1.)) * (1 - phi)
+            else:
+                return pseudo_equilibrium(
+                    K, phi, mol, T, n, gamma.f, gamma.args, 
+                    self.pseudo_equilibrium_inner_loop_options,
+                    self.pseudo_equilibrium_outer_loop_options,
+                )
         index = indices[-1]
         args = (mol, T, gamma.f, gamma.args)
         bounds = np.zeros([n, 2])

thermosteam/equilibrium/vle.py

@@ -1241,14 +1241,9 @@ def _refresh_v(self, V, y_bubble, x_dew, dz_bubble=None, dz_dew=None): # TODO: u
         self._v = v = self._estimate_v(V, y_bubble, x_dew, dz_bubble, dz_dew)
         self._V = V
         self._y = fn.normalize(v, v.sum() + self._F_mol_light)
-        try:
-            reload_cache = self._x is None or np.abs(self._z_last - self._z).sum() > 0.001
-        except:
-            reload_cache = True
-        if reload_cache:
-            l = self._mol_vle - v
-            l[l < 0] = 1e-12
-            self._x = fn.normalize(l, l.sum() + self._F_mol_heavy)
+        l = self._mol_vle - v
+        l[l < 0] = 1e-12
+        self._x = fn.normalize(l, l.sum() + self._F_mol_heavy)
 
     def _H_hat_err_at_T(self, T, H_hat, gas_conversion, liquid_conversion):
         v = self._solve_v(T, self._P, gas_conversion, liquid_conversion)

thermosteam/network.py

@@ -1426,6 +1426,30 @@ def add_bounded_numerical_specification(self, f=None, *args, **kwargs):
         self._specifications.append(BoundedNumericalSpecification(f, *args, **kwargs))
         return f
 
+    def run_phenomena(self):
+        """
+        Run mass and energy balance without converging phenomena. This method also runs specifications
+        user defined specifications unless it is being run within a 
+        specification (to avoid infinite loops). 
+        
+        See Also
+        --------
+        _run
+        specifications
+        add_specification
+        add_bounded_numerical_specification
+        
+        """
+        if self._skip_simulation_when_inlets_are_empty and all([i.isempty() for i in self._ins]):
+            for i in self._outs: i.empty()
+            return
+        if hasattr(self, '_run_phenomena'): 
+            self._run = self._run_phenomena
+            try: self._run_with_specifications()
+            finally: del self._run
+        else:
+            self._run_with_specifications()
+
     def run(self):
         """
         Run mass and energy balance. This method also runs specifications
@@ -1443,6 +1467,9 @@ def run(self):
         if self._skip_simulation_when_inlets_are_empty and all([i.isempty() for i in self._ins]):
             for i in self._outs: i.empty()
             return
+        self._run_with_specifications()
+
+    def _run_with_specifications(self):
         specifications = self._specifications
         if specifications:
             active_specifications = self._active_specifications
@@ -2066,9 +2093,12 @@ def fill_path(feed, path, paths_with_recycle,
     if not unit or unit._universal or unit not in units:
         paths_without_recycle.append(path)
     elif unit in path: 
-        path_with_recycle = path, feed
-        paths_with_recycle.append(path_with_recycle)
-        ends.add(feed)
+        if len(unit.outs) == 1 and unit.outs[0] in ends: 
+            paths_without_recycle.append(path)
+        else:
+            ends.add(feed)
+            path_with_recycle = path, feed
+            paths_with_recycle.append(path_with_recycle)
     elif feed in ends:
         paths_without_recycle.append(path)
     else:
@@ -2221,6 +2251,34 @@ def get_all_recycles(self, all_recycles=None):
             if isinstance(i, Network): i.get_all_recycles(all_recycles)
         return all_recycles
 
+    def sort_without_recycle(self, ends):
+        if self.recycle: return
+        path = tuple(self.path)
+
+        def sort(network):
+            subpath = tuple(network.path)
+            for o, j in enumerate(subpath):
+                if isinstance(j, Network):
+                    added = sort(j)
+                    if added: return True
+                elif j in sources: 
+                    self.path.remove(u)
+                    network.path.insert(o + 1, u)
+                    return True
+
+        for n, u in enumerate(path):
+            if isinstance(u, Network): continue
+            sources = set([i.source for i in u.ins if i.source and i not in ends])
+            subpath = path[n+1:]
+            for m, i in enumerate(subpath):
+                if isinstance(i, Network):
+                    added = sort(i)
+                    if added: break
+                elif i in sources: 
+                    self.path.remove(u)
+                    self.path.insert(n + m, u)
+                    break
+
     def sort(self, ends):
         isa = isinstance
         for i in self.path: 
@@ -2260,7 +2318,7 @@ def sort(self, ends):
         if not stop: warn(RuntimeWarning('network path could not be determined'))
 
     @classmethod
-    def from_feedstock(cls, feedstock, feeds=(), ends=None, units=None, final=True):
+    def from_feedstock(cls, feedstock, feeds=(), ends=None, units=None, final=True, recycles=True):
         """
         Create a Network object from a feedstock.
         
@@ -2290,9 +2348,10 @@ def from_feedstock(cls, feedstock, feeds=(), ends=None, units=None, final=True):
                 network.join_linear_network(linear_network) 
         else:
             network = Network([])
-        recycle_networks = [Network(*i) for i in cyclic_paths_with_recycle]
-        for recycle_network in recycle_networks:
-            network.join_recycle_network(recycle_network)
+        if recycles: 
+            recycle_networks = [Network(*i) for i in cyclic_paths_with_recycle]
+            for recycle_network in recycle_networks:
+                network.join_recycle_network(recycle_network)
         ends.update(network.streams)
         disjunction_streams = set([i.get_stream() for i in disjunctions])
         for feed in feeds:
@@ -2320,13 +2379,14 @@ def from_feedstock(cls, feedstock, feeds=(), ends=None, units=None, final=True):
             recycle_ends.update(disjunction_streams)
             recycle_ends.update(network.get_all_recycles())
             recycle_ends.update(tmo.utils.products_from_units(network.units))
+            network.sort_without_recycle(recycle_ends)
+            if recycles: network.reduce_recycles()
             network.sort(recycle_ends)
-            network.reduce_recycles()
             network.add_interaction_units()
         return network
 
     @classmethod
-    def from_units(cls, units, ends=None):
+    def from_units(cls, units, ends=None, recycles=True):
         """
         Create a System object from all units given.
 
@@ -2354,7 +2414,7 @@ def from_units(cls, units, ends=None):
             if not ends:
                 ends = tmo.utils.products_from_units(units) + [i.get_stream() for i in disjunctions]
             system = cls.from_feedstock(
-                feedstock, feeds, ends, units,
+                feedstock, feeds, ends, units, final=True, recycles=recycles,
             )
         else:
             system = cls(())