More Initializer Testing

idaes/core/initialization/general_hierarchical.py

@@ -66,14 +66,7 @@ def initialization_routine(
         Returns:
             Pyomo solver results object
         """
-        # Check to make sure model has something named "control_volume"
-        if not hasattr(model, "control_volume"):
-            raise TypeError(
-                f"Model {model.name} does not appear to be a standard form unit model. "
-                f"Please use an Initializer specific to the model being initialized."
-            )
-
-        # From here, the default initialization_routine is sufficient
+        # The default initialization_routine is sufficient
         return super().initialization_routine(
             model=model,
             plugin_initializer_args=plugin_initializer_args,
@@ -94,6 +87,13 @@ def initialize_main_model(self, model: Block, copy_inlet_state: bool = False):
             Pyomo solver results object from solve of main model
 
         """
+        # Check to make sure model has something named "control_volume"
+        if not hasattr(model, "control_volume"):
+            raise TypeError(
+                f"Model {model.name} does not appear to be a standard form unit model. "
+                f"Please use an Initializer specific to the model being initialized."
+            )
+
         # Get logger
         _log = self.get_logger(model)
 

idaes/core/initialization/initializer_base.py

@@ -151,6 +151,7 @@ def initialize(
         initial_guesses: dict = None,
         json_file: str = None,
         output_level=None,
+        exclude_unused_vars: bool = False,
     ):
         """
         Execute full initialization routine.
@@ -160,6 +161,7 @@ def initialize(
             initial_guesses: dict of initial guesses to load.
             json_file: file name of json file to load initial guesses from as str.
             output_level: (optional) output level to use during initialization run (overrides global setting).
+            exclude_unused_vars: whether to ignore unused variables when doing post-initialization checks.
 
         Note - can only provide one of initial_guesses or json_file.
 
@@ -196,7 +198,9 @@ def initialize(
             self._local_logger_level = None
 
         # 7. Check convergence
-        return self.postcheck(model, results_obj=results)
+        return self.postcheck(
+            model, results_obj=results, exclude_unused_vars=exclude_unused_vars
+        )
 
     def get_current_state(self, model: Block):
         """
@@ -638,7 +642,7 @@ def initialization_routine(
         self, model: Block, plugin_initializer_args: dict = None, **kwargs
     ):
         """
-        Common initialization routine for models with one control volume.
+        Common initialization routine for models with plugins.
 
         Args:
             model: Pyomo Block to be initialized
@@ -649,11 +653,6 @@ def initialization_routine(
         Returns:
             Pyomo solver results object
         """
-        if not hasattr(model, "control_volume"):
-            raise TypeError(
-                f"Model {model.name} does not appear to be a standard form unit model. "
-                f"Please use an Initializer specific to the model being initialized."
-            )
         if plugin_initializer_args is None:
             plugin_initializer_args = {}
 

idaes/core/initialization/tests/test_block_triangularization.py

@@ -16,11 +16,18 @@
 import pytest
 import types
 
-from pyomo.environ import ConcreteModel, Constraint, value, Var
+from pyomo.environ import ConcreteModel, Constraint, units, value, Var
+
+from idaes.core import FlowsheetBlock
 from idaes.core.initialization.block_triangularization import (
     BlockTriangularizationInitializer,
 )
 from idaes.core.initialization.initializer_base import InitializationStatus
+from idaes.models.unit_models.pressure_changer import (
+    Turbine,
+)
+
+from idaes.models.properties import iapws95
 
 __author__ = "Andrew Lee"
 

idaes/core/initialization/tests/test_general_hierarchical.py

@@ -32,6 +32,7 @@
 @pytest.mark.unit
 def test_nonstandard_model():
     m = ConcreteModel()
+    m.initialization_order = [m]
 
     initializer = SingleControlVolumeUnitInitializer()
 
@@ -116,7 +117,6 @@ def test_init_props_1D():
     m.control_volume.properties = Block()
 
     initializer = SingleControlVolumeUnitInitializer()
-    log = initializer.get_logger(m)
 
     def estimate_states(block, *args, **kwargs):
         block._estimated = True
@@ -240,9 +240,6 @@ def test_initialize_control_volume_copy():
     assert m.control_volume.reactions._initialized
 
 
-# TODO: init props 1D
-
-
 @pytest.mark.unit
 def test_initialize_submodels_no_order():
     m = ConcreteModel()

idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py

@@ -33,9 +33,10 @@
 
 References:
 
-1. "The properties of gases and liquids by Robert C. Reid"
-2. "Perry's Chemical Engineers Handbook by Robert H. Perry".
-3. H. Renon and J.M. Prausnitz, "Local compositions in thermodynamic excess functions for liquid mixtures.", AIChE Journal Vol. 14, No.1, 1968.
+#. "The properties of gases and liquids by Robert C. Reid"
+#. "Perry's Chemical Engineers Handbook by Robert H. Perry".
+#. H. Renon and J.M. Prausnitz, "Local compositions in thermodynamic excess
+   functions for liquid mixtures.", AIChE Journal Vol. 14, No.1, 1968.
 
 """
 # TODO: Missing docstrings
@@ -46,6 +47,7 @@
 
 # Import Pyomo libraries
 from pyomo.environ import (
+    Block,
     check_optimal_termination,
     Constraint,
     log,
@@ -58,7 +60,7 @@
     sqrt,
     units as pyunits,
 )
-from pyomo.common.config import ConfigValue, In
+from pyomo.common.config import ConfigDict, ConfigValue, In
 
 # Import IDAES cores
 from idaes.core import (
@@ -83,6 +85,7 @@
 from idaes.core.solvers import get_solver
 import idaes.core.util.scaling as iscale
 import idaes.logger as idaeslog
+from idaes.core.initialization import InitializerBase
 
 
 # Some more information about this module
@@ -301,12 +304,186 @@ def define_metadata(cls, obj):
         )
 
 
+class ActivityCoeffInitializer(InitializerBase):
+    """
+    Initializer for Activity Coefficient property packages.
+
+    This Initializer uses a hierarchical routine to initialize the
+    property package using the Pyomo solve_strongly_connected_components
+    method is used at each step to converge the problem.
+
+    Note that for systems without vapor-liquid equilibrium the generic
+    BlockTriangularizationInitializer is probably sufficient for initializing
+    the property package.
+
+    """
+
+    CONFIG = InitializerBase.CONFIG()
+    CONFIG.declare(
+        "solver",
+        ConfigValue(
+            default=None,
+            description="Solver to use for initialization",
+        ),
+    )
+    CONFIG.declare(
+        "solver_options",
+        ConfigDict(
+            implicit=True,
+            description="Dict of options to pass to solver",
+        ),
+    )
+    CONFIG.declare(
+        "calculate_variable_options",
+        ConfigDict(
+            implicit=True,
+            description="Dict of options to pass to 1x1 block solver",
+            doc="Dict of options to pass to calc_var_kwds argument in "
+            "solve_strongly_connected_components method. NOTE: models "
+            "involving ExternalFunctions must set "
+            "'diff_mode=differentiate.Modes.reverse_numeric'",
+        ),
+    )
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self._solver = None
+
+    def initialization_routine(
+        self,
+        model: Block,
+    ):
+        """
+        Sequential initialization routine for cubic EoS properties.
+
+        Args:
+            model: model to be initialized
+
+        Returns:
+            None
+        """
+        # Deactivate the constraints specific for outlet block i.e.
+        # when defined state is False
+        init_log = idaeslog.getInitLogger(
+            model.name, self.config.output_level, tag="properties"
+        )
+        solve_log = idaeslog.getSolveLogger(
+            model.name, self.config.output_level, tag="properties"
+        )
+
+        for k in model.values():
+            if (k.config.defined_state is False) and (
+                k.params.config.state_vars == "FTPz"
+            ):
+                k.eq_mol_frac_out.deactivate()
+
+        # Create solver
+        solver = get_solver(self.config.solver, self.config.solver_options)
+
+        # ---------------------------------------------------------------------
+        # Initialization sequence: Deactivating certain constraints
+        # for 1st solve
+        for k in model.values():
+            for c in k.component_objects(Constraint):
+                if c.local_name in [
+                    "eq_total",
+                    "eq_comp",
+                    "eq_phase_equilibrium",
+                    "eq_enth_mol_phase",
+                    "eq_entr_mol_phase",
+                    "eq_Gij_coeff",
+                    "eq_A",
+                    "eq_B",
+                    "eq_activity_coeff",
+                ]:
+                    c.deactivate()
+
+        # First solve for the active constraints that remain (p_sat, T_bubble,
+        # T_dew). Valid only for a 2 phase block. If single phase,
+        # no constraints are active.
+        # NOTE: "k" is the last value from the previous for loop
+        # only for the purpose of having a valid index. The assumption
+        # here is that for all values of "blk[k]", the attribute exists.
+        # TODO: This step can result in a non-square problem. This needs to be
+        # investigated
+        # TODO: Block Triangularization solver does not like this routine for some
+        # reason, maybe related to the above issue
+        if (
+            (k.config.has_phase_equilibrium)
+            or (k.config.parameters.config.valid_phase == ("Liq", "Vap"))
+            or (k.config.parameters.config.valid_phase == ("Vap", "Liq"))
+        ):
+
+            with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+                res = solve_indexed_blocks(solver, [model], tee=slc.tee)
+        else:
+            res = "skipped"
+        init_log.info(f"Initialization Step 1 {idaeslog.condition(res)}.")
+
+        # Continue initialization sequence and activate select constraints
+        for k in model.values():
+            for c in k.component_objects(Constraint):
+                if c.local_name in [
+                    "eq_total",
+                    "eq_comp",
+                    "eq_phase_equilibrium",
+                ]:
+                    c.activate()
+            if k.config.parameters.config.activity_coeff_model != "Ideal":
+                # assume ideal and solve
+                k.activity_coeff_comp.fix(1)
+
+        # Second solve for the active constraints
+        with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+            with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+                res = solve_indexed_blocks(solver, [model], tee=slc.tee)
+        init_log.info(f"Initialization Step 2 {idaeslog.condition(res)}.")
+
+        # Activate activity coefficient specific constraints
+        for k in model.values():
+            if k.config.parameters.config.activity_coeff_model != "Ideal":
+                for c in k.component_objects(Constraint):
+                    if c.local_name in ["eq_Gij_coeff", "eq_A", "eq_B"]:
+                        c.activate()
+
+        with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+            with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+                res = solve_indexed_blocks(solver, [model], tee=slc.tee)
+        init_log.info(f"Initialization Step 3 {idaeslog.condition(res)}.")
+
+        for k in model.values():
+            if k.config.parameters.config.activity_coeff_model != "Ideal":
+                k.eq_activity_coeff.activate()
+                k.activity_coeff_comp.unfix()
+
+        with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+            with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+                res = solve_indexed_blocks(solver, [model], tee=slc.tee)
+        init_log.info(f"Initialization Step 4 {idaeslog.condition(res)}.")
+
+        for k in model.values():
+            for c in k.component_objects(Constraint):
+                if c.local_name in ["eq_enth_mol_phase", "eq_entr_mol_phase"]:
+                    c.activate()
+
+        with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+            with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
+                res = solve_indexed_blocks(solver, [model], tee=slc.tee)
+        init_log.info(f"Initialization Step 5 {idaeslog.condition(res)}.")
+
+        return res
+
+
 class _ActivityCoeffStateBlock(StateBlock):
     """
     This Class contains methods which should be applied to Property Blocks as a
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    # Set default initializer
+    default_initializer = ActivityCoeffInitializer
+
     def fix_initialization_states(self):
         """
         Fixes state variables for state blocks.