implement tests for phenomenodes

Author: yoelcortes

biosteam/_system.py

@@ -2517,7 +2517,8 @@ def track_convergence(self, track=True):
         if track:
             self._track_convergence = True
             for i in self.stages: i.create_equation_nodes()
-            for i in self.stages: i.initialize_equation_nodes()
+            with self.stage_configuration(aggregated=False):
+                for i in self.stages: i.initialize_equation_nodes()
             equations = tuple(set(sum([i.equation_nodes for i in self.stages], ())))
             variables = tuple(set(sum([i.variables for i in equations], ())))
             self.variable_profiles = {i: [] for i in variables}
@@ -2641,7 +2642,9 @@ def _collect_subgraph_variables(self, configuration, subgraph):
                 f'only {all_subgraphs!r} are valid'
             )
         for i in configuration.nodes:
-            for variable, value in getattr(i, f'_collect_{subgraph}_variables')():
+            method = getattr(i, f'_collect_{subgraph}_variables', None)
+            if method is None: continue
+            for variable, value in method():
                 if variable in collected_variables: continue
                 variables[variable].append(value)
                 collected_variables.add(variable)

biosteam/units/stage.py

@@ -79,7 +79,11 @@ class SinglePhaseStage(Unit):
     _N_ins = 2
     _N_outs = 1
     _ins_size_is_fixed = False
-    _energy_variable = 'T'
+    
+    @property
+    def _energy_variable(self):
+        if self.T is None: return 'T'
+        else: return None
 
     def _init(self, T=None, P=None, Q=None, phase=None):
         self.T = T
@@ -156,7 +160,75 @@ def _update_energy_variable(self, departure):
         for i in self.outs: i.T += departure
 
     def _update_nonlinearities(self): pass
+    
+    @property
+    def equation_node_names(self): 
+        if self._energy_variable is None:
+            return (
+                'overall_material_balance_node', 
+            )
+        else:
+            return (
+                'overall_material_balance_node', 
+                'energy_balance_node',
+            )
+    
+    def initialize_overall_material_balance_node(self):
+        self.overall_material_balance_node.set_equations(
+            inputs=[j for i in self.ins if (j:=i.F_node)],
+            outputs=[i.F_node for i in self.outs],
+        )
+        
+    def initialize_energy_balance_node(self):
+        self.energy_balance_node.set_equations(
+            inputs=(
+                self.T_node, *[i.F_node for i in self.outs],
+                *[j for i in self.ins if (j:=i.E_node)],
+            ),
+            outputs=[
+                j for i in self.outs if (j:=i.E_node)
+            ],
+        )
+    
+    @property
+    def T_node(self):
+        if hasattr(self, '_T_node'): 
+            self._T_node.value = self.T
+            return self._T_node
+        self._T_node = var = VariableNode(f"{self.node_tag}.T",
+            self.T,
+            self.energy_balance_node,
+            *[i.sink.energy_balance_node for i in self.outs if i.sink],
+        )
+        return var 
+        
+    @property
+    def E_node(self):
+        if self._energy_variable is None:
+            return None
+        else:
+            return self.T_node
 
+    def get_connected_material_nodes(self, stream):
+        # TODO: Make this work for inlet streams and split streams
+        eqs = [self.overall_material_balance_node]
+        if self._energy_variable is not None:
+            eqs.append(self.energy_balance_node)
+        return eqs
+    
+    def _collect_material_balance_variables(self):
+        nodes = [i.F_node for i in self.outs]
+        return [
+            (i, i.value) for i in nodes
+        ]
+        # list[tuple[VariableNode, value]] of stage-name and material balance variable pairs
+    
+    def _collect_energy_balance_variables(self):
+        nodes = [j for i in self.outs if (j:=i.E_node)]
+        return [
+            (i, i.value) for i in nodes
+        ] # list[tuple[VariableNode, value]] of stage-name and energy balance variable pairs
+    
 
 class ReactivePhaseStage(bst.Unit): # Does not include VLE
     _N_outs = _N_ins = 1
@@ -253,7 +325,7 @@ def _update_nonlinearities(self):
         old = self.dmol
         new = self.reaction.conversion(self.ins[0])
         self.dmol = f * old + (1 - f) * new
-    
+
 
 # %% Two phases
 
@@ -479,7 +551,7 @@ def _run(self):
         self.partition._run()
         for i in self.splitters: i._run()
         self._update_separation_factors()
-        
+    
     def initialize_overall_material_balance_node(self):
         self.overall_material_balance_node.set_equations(
             inputs=[j for i in self.ins if (j:=i.F_node)],
@@ -561,7 +633,7 @@ def get_connected_material_nodes(self, stream):
                self.overall_material_balance_node]
         if stream.phase == 'l':
             eqs.append(self.phenomena_node)
-        elif self._energy_variable is None and stream.phase == 'g':
+        elif self._energy_variable is not None and stream.phase == 'g':
             eqs.append(self.energy_balance_node)
         return eqs
 

tests/test_phenomenode_tracking.py

@@ -10,20 +10,54 @@ def test_flash():
     import biosteam as bst
     import numpy as np
     from numpy.testing import assert_allclose
-    bst.System.strict_phenomena_decomposition = True
-    
+    bst.settings.set_thermo(['Water', 'Ethanol'], ideal=True)
+    feed_liquid = bst.Stream('feed_liquid', Water=50, Ethanol=50, phase='l', T = 400)
+    feed_gas = bst.Stream('feed_gas', Water=50, Ethanol=50, phase='g', T=400)
+    flash_1 = bst.StageEquilibrium('flash_1', ins=[feed_liquid, feed_gas], Q=0, P=101325, phases=('g', 'l'))
+    flashes = [flash_1]
+    for flash in flashes:
+        flash.T = 350
+        flash.K = np.array([0.5, 1.5])
+        flash.B = 1.
+        flash.outs[0].mol[:] = [50, 50]
+        flash.outs[1].mol[:] = [50, 50]
+    sys = bst.System('sys', path=flashes, algorithm='phenomena-oriented')
+    sys.track_convergence()
+    for i in range(10):
+        sys.run_phenomena()
+    phenomena_graph = sys.get_phenomena_graph()
+    shape = phenomena_graph.variable_profiles.shape
+    assert shape == (30, 8)
+        
 
 def test_2_stage_flash():
     import biosteam as bst
     import numpy as np
     from numpy.testing import assert_allclose
-    bst.System.strict_phenomena_decomposition = True
+    bst.settings.set_thermo(['Water', 'Ethanol'], ideal=True)
+    feed_liquid = bst.Stream('feed_liquid', Water=50, Ethanol=50, phase='l', T = 400)
+    feed_gas = bst.Stream('feed_gas', Water=50, Ethanol=50, phase='g', T=400)
+    flash_1 = bst.StageEquilibrium('flash_1', ins=[feed_liquid, feed_gas], Q=0, P=101325, phases=('g', 'l'))
+    flash_2 = bst.StageEquilibrium('flash_2', ins=[flash_1-0, feed_liquid], B=1, P=101325, phases=('g', 'l'))
+    flashes = [flash_1, flash_2]
+    for flash in flashes:
+        flash.T = 350
+        flash.K = np.array([0.5, 1.5])
+        flash.B = 1.
+        flash.outs[0].mol[:] = [50, 50]
+        flash.outs[1].mol[:] = [50, 50]
+    sys = bst.System('sys', path=flashes, algorithm='phenomena-oriented')
+    sys.track_convergence()
+    for i in range(10):
+        sys.run_phenomena()
+    phenomena_graph = sys.get_phenomena_graph()
+    shape = phenomena_graph.variable_profiles.shape
+    assert shape == (30, 15)
 
 def test_column():
     import biosteam as bst
     import numpy as np
     from numpy.testing import assert_allclose
-    bst.System.strict_phenomena_decomposition = True
     bst.settings.set_thermo(['Water', 'Ethanol'], cache=True)
     feed = bst.Stream('feed', Ethanol=80, Water=100, T=80.215 + 273.15)
     D1 = bst.MESHDistillation('D1', N_stages=5, ins=[feed], feed_stages=[2],
@@ -33,9 +67,88 @@ def test_column():
     D1.set_flow_rates(D1.hot_start())
     sys = bst.System('sys', path=[D1])
     sys.track_convergence(True)
-    for i in range(80):
+    for i in range(10):
+        sys.run_phenomena()
+    phenomena_graph = sys.get_phenomena_graph()
+    shape = phenomena_graph.variable_profiles.shape
+    assert shape == (30, 38)
+ 
+def test_system():
+    import biosteam as bst
+    import numpy as np
+    from numpy.testing import assert_allclose
+    bst.settings.set_thermo(['Water', 'AceticAcid', 'EthylAcetate'], cache=True)
+    solvent_feed_ratio = 1
+    @bst.SystemFactory
+    def system(ins, outs):
+        feed = bst.Stream('feed', AceticAcid=6660, Water=43600)
+        solvent = bst.Stream('solvent', EthylAcetate=65000)
+        recycle = bst.Stream('recycle')
+        LE = bst.MultiStageEquilibrium(
+            N_stages=6, ins=[feed, solvent, recycle],
+            feed_stages=(0, -1, -1),
+            phases=('L', 'l'),
+            maxiter=200,
+            use_cache=True,
+        )
+        # DAA = bst.MultiStageEquilibrium(N_stages=6, ins=[LE-0], feed_stages=[3],
+        #     outs=['vapor', 'liquid'],
+        #     stage_specifications={0: ('Reflux', 0.673), -1: ('Boilup', 2.57)},
+        #     maxiter=200,
+        #     phases=('g', 'l'),
+        #     use_cache=True,
+        # )
+        DEA = bst.MultiStageEquilibrium(N_stages=6, ins=[LE-1], feed_stages=[3],
+            outs=['vapor', 'liquid'],
+            stage_specifications={0: ('Reflux', 0.673), -1: ('Boilup', 2.57)},
+            phases=('g', 'l'),
+            maxiter=200,
+            use_cache=True,
+        )
+        HX = bst.SinglePhaseStage(ins=DEA-0, outs=recycle, T=320, phase='l')
+        
+        chemicals = bst.settings.chemicals
+        
+        @LE.add_specification(run=True)
+        def fresh_solvent_flow_rate():
+            broth = feed.F_mass
+            EtAc_recycle = recycle.imass['EthylAcetate']
+            EtAc_required = broth * solvent_feed_ratio
+            if EtAc_required < EtAc_recycle:
+                recycle.F_mass *= EtAc_required / EtAc_recycle
+                EtAc_recycle = recycle.imass['EthylAcetate']
+            EtAc_fresh = EtAc_required - EtAc_recycle
+            solvent.imass['EthylAcetate'] = max(
+                0, EtAc_fresh
+            )
+        
+        @solvent.material_balance
+        def fresh_solvent_flow_rate():
+            s = np.ones(chemicals.size)
+            r = np.zeros(chemicals.size)
+            v = r.copy()
+            index = chemicals.index('EthylAcetate')
+            r[index] = 1
+            v[index] = solvent_feed_ratio * feed.F_mass / chemicals.EthylAcetate.MW
+            return (
+                {solvent: s,
+                 recycle: r},
+                 v
+            )
+    sys = system(algorithm='phenomena oriented', 
+                    molar_tolerance=1e-9,
+                    relative_molar_tolerance=1e-9,
+                    maxiter=2000,
+                    method='fixed-point')
+    sys.track_convergence(True)
+    for i in range(10):
         sys.run_phenomena()
     phenomena_graph = sys.get_phenomena_graph()
+    shape = phenomena_graph.variable_profiles.shape
+    assert shape == (37, 133)
 
 if __name__ == '__main__':
-    pass
+    test_flash()
+    test_2_stage_flash()
+    test_column()
+    test_system()