solve merge conflict

Author: yoelcortes

biosteam/_settings.py

@@ -165,22 +165,23 @@ def define_allocation_property(
         self, name: str, basis: float, 
         stream: Callable=None, 
         power_utility: Callable=None,
-        heat_utility: Callable=None
+        heat_utility: Callable=None,
+        safe: bool=True,
     ):
     """Define a new allocation property by property getters."""
     allocation_name = name + '-allocation'
     units = basis + '/hr'
     if stream is not None:
         bst.Stream.define_property(
-            allocation_name, units, stream,
+            allocation_name, units, stream, safe=safe
         )
     if power_utility is not None:
         bst.PowerUtility.define_property(
-            allocation_name, units, power_utility,
+            allocation_name, units, power_utility, safe=safe
         )
     if heat_utility is not None:
         bst.HeatUtility.define_property(
-            allocation_name, units, heat_utility,
+            allocation_name, units, heat_utility, safe=safe
         )
     bst.allocation_properties[name] = basis
 
@@ -255,16 +256,19 @@ def revenue_allocation(stream):
     stream=lambda self: max(self.LHV, 0),
     power_utility=lambda self: max(-self.rate * 3600, 0.),
     heat_utility=lambda self: max(self.duty, 0),
+    safe=False,
 )
 settings.define_allocation_property(
     'revenue', 'USD', 
     stream=revenue_allocation,
     power_utility=lambda self: max(-self.cost, 0.),
     heat_utility=lambda self: max(self.cost, 0),
+    safe=False,
 )
 settings.define_allocation_property(
-    'mass', 'kg', stream=lambda self: np.dot(self.chemicals.MW, self.mol)
+    'mass', 'kg', stream=lambda self: np.dot(self.chemicals.MW, self.mol),
+    safe=False,
 )
 settings.define_impact_indicator(
-    'WU', 'L' # Water usage is a special indicator native to BioSTEAM
+    'WU', 'L', # Water usage is a special indicator native to BioSTEAM
 )

biosteam/_system.py

@@ -2682,9 +2682,8 @@ def track_convergence(
                         'energy': [],
                         'phenomenode': [],
                     }
-            if not self.recycle:
-                self.timer = Timer()
-                self.timer.start()
+            self.timer = Timer()
+            self.timer.start()
 
         else:
             for i in self.variable_profiles: i.getter = None
@@ -2696,6 +2695,7 @@ def track_convergence(
             self.edge_profiles = None
             self.grouped_variables = None
             for i in self.units: i.tracking = False
+            self.timer = None
 
     def get_profiles(self, equation_profiles, edge_profiles):
         time, variable_profiles = self.get_variable_profiles()
@@ -2738,9 +2738,9 @@ def get_phenomegraph(self, bipartite=False, **kwargs):
     def get_monopartite_phenomegraph(self, decomposition=None, dotfile=None):
         if decomposition is None: decomposition = self.algorithm.lower()
         subgraph_units = [i.ID for i in self.path]
-        if 'phenomenode' in decomposition:
+        if 'phenomena' in decomposition:
             if 'modular' in decomposition:
-                criteria = [*all_subgraphs, *[(i, 'phenomenode') for i in subgraph_units]]
+                criteria = [*all_subgraphs, *[(i, 'phenomena') for i in subgraph_units]]
             else:
                 criteria = all_subgraphs
         elif 'modular' in decomposition:
@@ -2804,7 +2804,7 @@ def get_bipartite_phenomegraph(
             raise ValueError('unknown decompostion')
 
         def equation_match(equation, subgraph):
-            name = equation.name.replace('phenomenode', 'phenomena')
+            name = equation.name
             if isinstance(subgraph, tuple):
                 return all([i in name for i in subgraph])
             else:

biosteam/process_tools/system_factory.py

@@ -237,9 +237,6 @@ def __call__(self, ID=None, ins=None, outs=None,
                 bst.settings.set_thermo(fthermo(chemicals=chemicals))
             elif not hasattr(bst.settings, '_thermo'):
                 bst.settings.set_thermo(fthermo())
-        if autorename is not None: 
-            original_autorename = tmo.utils.Registry.AUTORENAME
-            tmo.utils.Registry.AUTORENAME = autorename
         ins = create_streams(self.ins, ins, 'inlets', self.fixed_ins_size)
         outs = create_streams(self.outs, outs, 'outlets', self.fixed_outs_size)
         rename = area is not None
@@ -275,7 +272,13 @@ def __call__(self, ID=None, ins=None, outs=None,
             ins = tuple([module.auxin(i) for i in module.ins])
             outs = tuple([module.auxout(i) for i in module.outs])
             with bst.Flowsheet(ID), bst.System(**options) as system:
-                self.f(ins, outs, **kwargs)
+                if autorename is None:
+                    self.f(ins, outs, **kwargs)
+                else:
+                    original_autorename = tmo.utils.Registry.AUTORENAME
+                    tmo.utils.Registry.AUTORENAME = autorename
+                    try: self.f(ins, outs, **kwargs)
+                    finally: tmo.utils.Registry.AUTORENAME = original_autorename
             module._init(system=system)
         else: # Create system
             with (bst.MockSystem() if mockup else bst.System(**options)) as system:
@@ -286,10 +289,15 @@ def __call__(self, ID=None, ins=None, outs=None,
                 elif udct:
                     unit_registry = system.flowsheet.unit
                     irrelevant_units = set(unit_registry)
-                self.f(ins, outs, **kwargs)
+                if autorename is None:
+                    self.f(ins, outs, **kwargs)
+                else:
+                    original_autorename = tmo.utils.Registry.AUTORENAME
+                    tmo.utils.Registry.AUTORENAME = autorename
+                    try: self.f(ins, outs, **kwargs)
+                    finally: tmo.utils.Registry.AUTORENAME = original_autorename
         system.load_inlet_ports(ins, {k: i for i, j in enumerate(self.ins) if (k:=get_name(j)) is not None})
         system.load_outlet_ports(outs, {k: i for i, j in enumerate(self.outs) if (k:=get_name(j)) is not None})
-        if autorename is not None: tmo.utils.Registry.AUTORENAME = original_autorename
         if udct: 
             unit_dct = {}
             def add(key, unit):

biosteam/units/abstract_stirred_tank_reactor.py

@@ -269,6 +269,9 @@ class AbstractStirredTankReactor(PressureVessel, Unit, isabstract=True):
     #: Phases accounted for in residence time
     tau_phases_default: str = 'l'
 
+    #: Basis of residence time. Defaults to inlet flow rates
+    tau_basis: str = 'ins' # Alternatively, 'outs'
+    
     @property
     def effluent(self):
         return self.outs[-1]
@@ -426,12 +429,14 @@ def _get_duty(self): # User can replace
 
     def _design(self, size_only=False):
         Design = self.design_results
-        ins_F_vol = sum([i.F_vol for i in self.ins if i.phase in self.tau_phases])
+        tau_basis = self.tau_basis
+        streams = getattr(self, tau_basis)
+        F_vol = sum([i.F_vol for i in streams if i.phase in self.tau_phases])
         P_pascal = (self.P if self.P else self.outs[0].P)
         P_psi = P_pascal * 0.000145038 # Pa to psi
         length_to_diameter = self.length_to_diameter
         if self.batch:
-            v_0 = ins_F_vol
+            v_0 = F_vol
             tau = self.tau
             tau_0 = self.tau_0
             V_wf = self.V_wf
@@ -446,14 +451,14 @@ def _design(self, size_only=False):
             N = Design['Number of reactors']
         else:
             if self.N is None:
-                V_total = ins_F_vol * self.tau / self.V_wf
+                V_total = F_vol * self.tau / self.V_wf
                 N = ceil(V_total/self.V_max)
                 if N == 0:
                     V_reactor = 0
                 else:
                     V_reactor = V_total / N
             else:
-                V_total = ins_F_vol * self.tau / self.V_wf
+                V_total = F_vol * self.tau / self.V_wf
                 V_reactor = V_total / self.N
             Design['Reactor volume'] = V_reactor
             Design['Number of reactors'] = N

biosteam/units/design_tools/Gibbs_equilibrium_reaction.py

@@ -16,37 +16,54 @@ def Gibbs_equilibrium_objective(
         product, # Product stream (constant T and P)
         main_phase, # Initial guess phase
         phase_equilibrium, # Name of phase equilibrium. Valid inputs include None, 'vle', 'lle', and 'vlle'
+        phase_hook,
     ):
     mask = mass < 0
     penalty = mask.any()
     if penalty: mass[mask] = 0
     product.empty()
     product[main_phase].imol[IDs] = mass / MWs
-    if phase_equilibrium is not None:
-        eq = getattr(product, phase_equilibrium)
-        eq(T=eq.T, P=eq.P)
+    if phase_hook is None:
+        if phase_equilibrium is not None:
+            eq = getattr(product, phase_equilibrium)
+            eq(T=eq.T, P=eq.P)
+    else:
+        phase_hook(product)
     return product.G
 
 def minimize_Gibbs_free_energy(
         product, 
         IDs=None, # Names of potential products
         method=None,
+        phase_hook=None,
     ):
     if method is None: method = 'differential evolution'
+    
+    # Normalize to 1 kg/hr
     F_mass = product.F_mass
-    mol_norm = product.mol / F_mass # Normalize
+    mol_norm = product.mol / F_mass
+    
+    # Get atomic flows 
     chemicals = product.chemicals
     atoms = chemicals.formula_array @ mol_norm 
     index, = np.where(atoms)
     atoms = atoms[index]
     formula_array = chemicals.formula_array[index]
+    
+    # Reduce to only possible products (such that all atoms exist in feed)
     if IDs is None: 
         IDs = chemicals.IDs
         MWs = chemicals.MW
     else:
         index = chemicals.get_index(IDs)
         formula_array = formula_array[:, index]
         MWs = chemicals.MW[index]
+    index, = np.where(formula_array.any(axis=0))
+    IDs = [IDs[i] for i in index]
+    formula_array = formula_array[:, index]
+    MWs = MWs[index]
+    
+    # Specify atomic and mass constraints
     N = len(IDs)
     lb = np.zeros(N)
     ub = np.ones(N)
@@ -59,24 +76,41 @@ def minimize_Gibbs_free_energy(
     atomic_balance = LinearConstraint(
         formula_array / MWs, atoms, atoms
     )
-    match product.phases:
-        case ('g', 'l'):
-            phase_equilibrium = 'vle' 
+    
+    # Choose main phase and phase equilibrium algorithm
+    if phase_hook:
+        phases = product.phases
+        if 'g' in phases:
             main_phase = 'g'
-        case ('L', 'l'):
-            phase_equilibrium = 'lle' 
-            main_phase = 'l'
-        case ('L', 'g', 'l'):
-            phase_equilibrium = 'vlle'
+        elif 'l' in phases:
             main_phase = 'l'
-        case ('l', 's'):
-            phase_equilibrium = 'sle'
-            main_phase = 'l'
-        case [main_phase]:
-            phase_equilibrium = None
-        case _:
-            raise RuntimeError(f'phase equilibrium for {product.phases!r} not supported')
-    f_args = (MWs, IDs, product, main_phase, phase_equilibrium)
+        elif 'L' in phases:
+            main_phase = 'L'
+        elif 's' in phases:
+            main_phase = 's'
+        else:
+            raise RuntimeError('main phase could not be found')
+    else:
+        match product.phases:
+            case ('g', 'l'):
+                phase_equilibrium = 'vle' 
+                main_phase = 'g'
+            case ('L', 'l'):
+                phase_equilibrium = 'lle' 
+                main_phase = 'l'
+            case ('L', 'g', 'l'):
+                phase_equilibrium = 'vlle'
+                main_phase = 'l'
+            case ('l', 's'):
+                phase_equilibrium = 'sle'
+                main_phase = 'l'
+            case [main_phase]:
+                phase_equilibrium = None
+            case _:
+                raise RuntimeError(f'phase equilibrium for {product.phases!r} not supported')
+                
+    # Solve
+    f_args = (MWs, IDs, product, main_phase, phase_equilibrium, phase_hook)
     if method == 'differential evolution':
         polish = lambda *args, **kwargs: minimize(
             *args, **kwargs, 
@@ -89,7 +123,7 @@ def minimize_Gibbs_free_energy(
             bounds=bounds,
             constraints=[atomic_balance],
             polish=polish,
-            tol=1e-5,
+            tol=1e-6,
             seed=0,
         )
     elif method == 'COBYLA':
@@ -107,9 +141,13 @@ def minimize_Gibbs_free_energy(
             "only 'SHGO' and 'differential evolution' "
             "are valid methods"
         )
+    
+    # Set solution
     product.empty()
     product[main_phase].imol[IDs] = solution.x / MWs
-    if phase_equilibrium is not None:
+    if phase_hook:
+        phase_hook(product)
+    elif phase_equilibrium is not None:
         eq = getattr(product, phase_equilibrium)
         eq(T=eq.T, P=eq.P)
     product.F_mass = F_mass

biosteam/units/equilibrium_reactor.py

@@ -8,7 +8,7 @@
 import biosteam as bst
 from .abstract_stirred_tank_reactor import AbstractStirredTankReactor
 from .design_tools.Gibbs_equilibrium_reaction import minimize_Gibbs_free_energy
-from typing import Optional, Iterable
+from typing import Optional, Iterable, Callable
 
 __all__ = (
     'EquilibriumReactor',
@@ -153,18 +153,21 @@ class EquilibriumReactor(AbstractStirredTankReactor):
     kW_per_m3_default = 0
     batch_default = False
     tau_phases_default = 'lg'
+    tau_basis = 'outs' # Alternatively, 'ins'
     method_default = 'differential evolution' # Alternatively 'COBYLA'
 
     def _init(self, 
             phases: Optional[str]=None, 
             products: Optional[Iterable[str]]=None, 
             method: Optional[str]=None,
+            phase_hook: Optional[Callable]=None,
             **kwargs
         ):
         AbstractStirredTankReactor._init(self, **kwargs)
         self.phases = phases
         self.products = products
         self.method = self.method_default if method is None else method
+        self.phase_hook = phase_hook
 
     def _get_duty(self): return self.Hnet
 
@@ -186,5 +189,5 @@ def _run(self):
         product.P = self.P
         product.empty()
         outs[0].mix_flows(ins)
-        minimize_Gibbs_free_energy(product, self.products, method=self.method)
+        minimize_Gibbs_free_energy(product, self.products, method=self.method, phase_hook=self.phase_hook)