2.52.8 release; improvements to phenomena-based sim

Author: yoelcortes

biosteam/__init__.py

@@ -13,7 +13,7 @@
 
 """
 from __future__ import annotations
-__version__ = '2.52.7'
+__version__ = '2.52.8'
 
 #: Chemical engineering plant cost index (defaults to 567.5 at 2017).
 CE: float = 567.5 

biosteam/_system.py

@@ -94,8 +94,6 @@ class Configuration:
                  'composition_sensitive_nodes',
                  '_has_dynamic_coefficients')
 
-    energy_variable_relaxation_factor = 0.5
-    
     def __init__(self, path, stages, nodes, streams, stream_ref, connections, aggregated):
         self.path = path
         self.stages = stages
@@ -139,24 +137,22 @@ def solve_nonlinearities(self):
                     i._update_nonlinearities()
         else:
             for i in self.path: 
-                if getattr(i, 'decoupled', False): 
-                    assert False
-                    i.run()
                 i._update_nonlinearities()
 
         if self.composition_sensitive_path:
             # for i in self.composition_sensitive_path: i._run()
+            # breakpoint()
             containers, flows = self.solve_material_flows(composition_sensitive=True, full_output=True)
             def update_inner_material_balance_parameters(flows):
-                for i in self.composition_sensitive_path: i._update_composition_parameters()
+                for i in self.composition_sensitive_nodes: i._update_composition_parameters()
                 return self.solve_material_flows(composition_sensitive=True)
 
             flx.fixed_point(
                 update_inner_material_balance_parameters,
                 flows, 
                 xtol=tmo.LLE.pseudo_equilibrium_inner_loop_options['xtol'], 
                 rtol=tmo.LLE.pseudo_equilibrium_inner_loop_options['rtol'], 
-                maxiter=10, checkiter=False,
+                maxiter=100, checkiter=False,
                 checkconvergence=False,
             )
             for i in self.composition_sensitive_path: i._update_net_flow_parameters()
@@ -243,17 +239,6 @@ def solve_material_flows(self, composition_sensitive=False, full_output=False):
                 values[values < 0] = 0
                 for mol, value in zip(containers, values): # update material flows
                     mol[:] = value
-                # for obj, mask, value in zip(objs, masks, values): # update material flows
-                #     indexer, phase = obj
-                #     index = indexer._phase_indexer(phase)
-                #     mol = indexer.data.rows[index]
-                #     if mask.any():
-                #         new_mol = value[mask]
-                #         relaxation_factor = (-new_mol / (mol[mask] - new_mol)).max()
-                #         print(relaxation_factor)
-                #         mol[:] = relaxation_factor * mol + (1 - relaxation_factor) * value
-                #     else:
-                #         mol[:] = value
             else:
                 if full_output:
                     containers = []
@@ -276,37 +261,52 @@ def solve_material_flows(self, composition_sensitive=False, full_output=False):
         else:
             return values
 
-    def solve_energy_departures(self):
+    def solve_energy_flows(self):
         nodes = self.nodes
         A = []
         b = []
         for node in nodes:
-            for coefficients, value in node._create_energy_departure_equations():
+            for coefficients, value in node._create_energy_balance_equations():
+                A.append(coefficients)
+                b.append(value)
+        if not A: return
+        for node in nodes:
+            for coefficients, value in node._create_bulk_balance_equations():
                 A.append(coefficients)
                 b.append(value)
+        A = [{(getattr(i, 'material_reference', i), v): j for (i, v), j in x.items()} for x in A]
         A, objs = dictionaries2array(A)
-        departures = solve(A, np.array(b).T).T
-        # Prevent negative numbers
-        f = 1
-        for obj, delta in zip(objs, departures): 
-            if delta >= 0: continue
-            var = obj._energy_variable
-            value = getattr(obj, var, None)
-            if value is None: 
-                if var == 'T':
-                    value = obj.outs[0].T
+        values = solve(A, np.array(b).T).T
+        # try: values = solve(A, np.array(b).T).T
+        # except: 
+        #     A = []
+        #     b = []
+        #     breakpoint()
+        #     for node in nodes:
+        #         for coefficients, value in node._create_energy_balance_equations():
+        #             breakpoint()
+        #             A.append(coefficients)
+        #             b.append(value)
+        #     if not A: return
+        #     for node in nodes:
+        #         for coefficients, value in node._create_bulk_balance_equations():
+        #             breakpoint()
+        #             A.append(coefficients)
+        #             b.append(value)
+        for (obj, var), value in zip(objs, values):
+            if var == 'F_mol':
+                indexer, phase = obj
+                index = indexer._phase_indexer(phase)
+                mol = indexer.data.rows[index]
+                if value == 0:
+                    mol[:] = 0
                 else:
-                    raise RuntimeError(f'energy variable {var!r} value is None')
-            f = min(- value / delta, f)
-        departures *= (1 - self.energy_variable_relaxation_factor) * f
-        try:
-            for obj, departure in zip(objs, departures):
-                obj._update_energy_variable(departure)
-        except AttributeError as e:
-            if obj._update_energy_variable:
-                raise e
+                    total = mol.sum()
+                    if total != 0: mol[:] *= value / total
             else:
-                raise NotImplementedError(f'{obj!r} has no method `_update_energy_variable`')
+                obj._update_variable(var, value)
+        for node in nodes:
+            if hasattr(node, '_reset_bulk_variable'): node._reset_bulk_variable()
 
     def __enter__(self):
         units = self.stages
@@ -2526,13 +2526,14 @@ def run_phenomena(self):
             try:
                 conf.solve_nonlinearities()
                 if flag: self._collect_variables('phenomenode')
-                conf.solve_energy_departures()
+                conf.solve_energy_flows()
                 if flag: self._collect_variables('energy')
                 conf.solve_material_flows()
                 if flag: self._collect_variables('material')
             except (NotImplementedError, UnboundLocalError, KeyError) as error:
                 raise error
             except Exception as e:
+                raise e
                 print('FAILED!')
                 print(e)
                 print('-------')
@@ -2553,7 +2554,7 @@ def run_phenomena_modular(self):
                     if flag: self._collect_variables((i.ID, 'phenomenode'))
                     conf.solve_material_flows()
                     if flag: self._collect_variables('material')
-                    conf.solve_energy_departures()
+                    conf.solve_energy_flows()
                     if flag: self._collect_variables('energy')
             except (NotImplementedError, UnboundLocalError, KeyError) as error:
                 raise error

biosteam/_unit.py

@@ -157,10 +157,11 @@ def __init_subclass__(cls,
         if default_phenomena is None:
             methods = (
                 '_update_nonlinearities',
-                '_get_energy_departure_coefficient',
+                '_update_energy_coefficient',
+                '_update_variable',
                 '_create_material_balance_equations',
-                '_create_energy_departure_equations',
-                '_update_energy_variable',
+                '_create_energy_balance_equations',
+                '_create_bulk_balance_equations',
                 '_energy_variable',
             )
             if all([getattr(cls, i) is getattr(Unit, i) for i in methods]):
@@ -377,11 +378,11 @@ def _update_nonlinearities(self):
         self._duty = self.Hnet
         Ts_new = [i.T for i in outs]
         if all([i == j for i, j in zip(Ts_new, Ts)]): # T constant
-            self._energy_variable = None
+            self.specified_variable = 'T'
         elif self.Hnet == Q: # Q constant    
-            self._energy_variable = 'T'
+            self.specified_variable = 'Q'
         else:
-            self._energy_variable = None
+            self.specified_variable = 'B'
         for i, j in zip(outs, data): i.set_data(j)
 
     def _simulation_error(self):
@@ -401,20 +402,35 @@ def _simulation_error(self):
         self._outs = outs
         return np.abs(new_flows - flows).sum(), np.abs(Ts_new - Ts).sum()
 
-    def _get_energy_departure_coefficient(self, stream):
+    def _update_energy_coefficient(self, stream, coefficients):
         """
-        tuple[object, float] Return energy departure coefficient of a stream 
-        for phenomena-based simulation.
+        Update coefficient of a stream 
+        for phenomena-based simulation and 
+        return the reference value.
         """
-        if self._energy_variable == 'T': return (self, stream.C)
+        if self.specified_variable == 'Q': 
+            C = stream.C
+            coefficients[self, 'T'] = -C
+            return -C * stream.T
+        return 0
+    
+    def _update_variable(self, variable, value):
+        """
+        Update variable being solved in equations for 
+        phenomena-based simulation.
+        """
+        if variable == 'T':
+            for i in self.outs: i.T = value
+        else:
+            raise ValueError('invalid variable')
 
     def _create_material_balance_equations(self, composition_sensitive):
         """
         list[tuple[dict, array]] Create material balance equations for 
         phenomena-based simulation.
         """
         if self._link_streams: return []
-        fresh_inlets, process_inlets, equations = self._begin_equations(composition_sensitive)
+        fresh_inlets, process_inlets, equations = self._begin_material_equations(composition_sensitive)
         outs = self.flat_outs
         N = self.chemicals.size
         ones = np.ones(N)
@@ -425,7 +441,7 @@ def _create_material_balance_equations(self, composition_sensitive):
         except:
             rhs = predetermined_flow
         mol_total = sum([i.mol for i in outs])
-        for i, s in enumerate(outs):
+        for s in outs:
             split = s.mol / mol_total
             minus_split = -split
             eq_outs = {}
@@ -436,31 +452,58 @@ def _create_material_balance_equations(self, composition_sensitive):
             )
         return equations
 
-    def _create_energy_departure_equations(self):
+    def _create_energy_balance_equations(self):
         """
         list[tuple[dict, float]] Create energy departure equations for 
         phenomena-based simulation.
         """
+        fresh_inlets, process_inlets, equations = self._begin_energy_equations()
         if self._energy_variable == 'T':
-            coeff = {self: sum([i.C for i in self.outs])}
-            for i in self.ins: i._update_energy_departure_coefficient(coeff)
-            if self._dmol.any():
-                dH = self._duty - self.Hnet
+            dmol = self._dmol
+            if dmol.any(): 
+                dH = self._duty + sum([i.Hf for i in self.outs]) - sum([i.Hf for i in self.ins]) 
             else:
-                dH = self._duty + self.H_in - self.H_out
-            return [(coeff, dH)]
+                dH = self._duty
+            coeff = {(self, 'T'): sum([i.C for i in self.outs])}
+            for i, s in enumerate(self.outs): coeff[self, i] = s.h
+            for i in process_inlets: 
+                dH += i._update_energy_coefficient(coeff)
+            for i in fresh_inlets:
+                dH += i.H
+            equations.append(
+                (coeff, dH)
+            )
+            return equations
         else:
-            return []
+            return equations
 
-    def _update_energy_variable(self, departure):
+    def _create_bulk_balance_equations(self):
         """
-        Update energy variable being solved in energy departure equations for 
+        list[tuple[dict, array]] Create bulk balance equations for 
         phenomena-based simulation.
         """
-        if self._energy_variable == 'T':
-            for i in self.outs: i.T += departure
+        if self._link_streams: return []
+        fresh_inlets, process_inlets, equations = self._begin_bulk_equations()
+        outs = self.flat_outs
+        predetermined_flow = [i.F_mol for i in fresh_inlets]
+        try:
+            dF_mol = self._dmol.sum()
+            rhs = predetermined_flow + dF_mol
+        except:
+            rhs = predetermined_flow
+        F_mol_total = sum([i.F_mol for i in outs])
+        for i, s in enumerate(outs):
+            split = s.F_mol / F_mol_total
+            minus_split = -split
+            eq_outs = {}
+            for i in process_inlets: eq_outs[i, 'F_mol'] = minus_split
+            eq_outs[s, 'F_mol'] = 1
+            equations.append(
+                (eq_outs, split * rhs)
+            )
+        return equations
 
-    def _begin_equations(self, composition_sensitive):
+    def _begin_material_equations(self, composition_sensitive):
         inlets = self.ins
         fresh_inlets = []
         process_inlets = []
@@ -519,6 +562,56 @@ def _begin_equations(self, composition_sensitive):
             equations = material_equations
         return fresh_inlets, process_inlets, equations
 
+    def _begin_bulk_equations(self):
+        inlets = self.ins
+        fresh_inlets = []
+        process_inlets = []
+        equations = [f() for i in inlets for f in i.bulk_equations]
+        isfeed = lambda s: s.isfeed() or s.source._recycle_system is not self._recycle_system
+        if equations:
+            for i in inlets: 
+                if (isfeed(i) and 
+                    not i.bulk_equations):
+                    fresh_inlets.append(i)
+                elif len(i) > 1:
+                    process_inlets.extend(i)
+                else:
+                    process_inlets.append(i)  
+        else:
+             for i in inlets: 
+                 if isfeed(i):
+                     fresh_inlets.append(i)
+                 elif len(i) > 1:
+                     process_inlets.extend(i)
+                 else:
+                     process_inlets.append(i)
+        return fresh_inlets, process_inlets, equations
+    
+    def _begin_energy_equations(self):
+        inlets = self.ins
+        fresh_inlets = []
+        process_inlets = []
+        equations = [f() for i in inlets for f in i.energy_equations]
+        isfeed = lambda s: s.isfeed() or s.source._recycle_system is not self._recycle_system
+        if equations:
+            for i in inlets: 
+                if (isfeed(i) and 
+                    not i.energy_equations):
+                    fresh_inlets.append(i)
+                elif len(i) > 1:
+                    process_inlets.extend(i)
+                else:
+                    process_inlets.append(i)  
+        else:
+             for i in inlets: 
+                 if isfeed(i):
+                     fresh_inlets.append(i)
+                 elif len(i) > 1:
+                     process_inlets.extend(i)
+                 else:
+                     process_inlets.append(i)
+        return fresh_inlets, process_inlets, equations    
+    
     Inlets = piping.Inlets
     Outlets = piping.Outlets