Enable RO0D to be dynamic

.github/workflows/checks.yml

@@ -119,7 +119,7 @@ jobs:
         pip show idaes-pse
         echo '::endgroup::'
         echo '::group::Output of "idaes get-extensions" command'
-        idaes get-extensions --verbose
+        idaes get-extensions --extra petsc --verbose
         echo '::endgroup::'
     - name: Add coverage report pytest options
       if: matrix.coverage
@@ -218,7 +218,7 @@ jobs:
         pip show idaes-pse
         echo '::endgroup::'
         echo '::group::Output of "idaes get-extensions" command'
-        idaes get-extensions --verbose
+        idaes get-extensions --extra petsc --verbose
         echo '::endgroup::'
     - name: Run pytest
       run: |
@@ -263,7 +263,7 @@ jobs:
         pip show idaes-pse
         echo '::endgroup::'
         echo '::group::Output of "idaes get-extensions" command'
-        idaes get-extensions --verbose
+        idaes get-extensions --extra petsc --verbose
         echo '::endgroup::'
     - name: Install Jupyter kernel
       run: |

watertap/core/membrane_channel_base.py

@@ -108,7 +108,7 @@ class FrictionFactor(Enum):
     "dynamic",
     ConfigValue(
         default=False,
-        domain=In([False]),
+        domain=Bool,
         description="Dynamic model flag - must be False",
         doc="""Indicates whether this model will be dynamic or not.
 **default** - False. Membrane units do not yet support dynamic
@@ -120,7 +120,7 @@ class FrictionFactor(Enum):
     "has_holdup",
     ConfigValue(
         default=False,
-        domain=In([False]),
+        domain=Bool,
         description="Holdup construction flag - must be False",
         doc="""Indicates whether holdup terms should be constructed or not.
 **default** - False. Membrane units do not have defined volume, thus
@@ -1011,7 +1011,17 @@ def calculate_scaling_factors(self):
 
 # helper for validating configuration arguments for this CV
 def validate_membrane_config_args(unit):
-
+    if unit.config.dynamic and unit.config.has_holdup:
+        if not (
+            (unit.config.pressure_change_type != PressureChangeType.fixed_per_stage)
+            or (
+                unit.config.mass_transfer_coefficient
+                == MassTransferCoefficient.calculated
+            )
+        ):
+            raise ConfigurationError(
+                f"dynamics=True and has_holdup=True will not work while pressure_change_type={unit.config.pressure_change_type} and mass_tranfer_coefficient={unit.config.mass_transfer_coefficient}\n. To enable dynamics, choose any other configuration option for pressure_change_type or mass_transfer_coefficient."
+            )
     if (
         unit.config.pressure_change_type is not PressureChangeType.fixed_per_stage
         and unit.config.has_pressure_change is False

watertap/property_models/NaCl_prop_pack.py

@@ -99,6 +99,14 @@ def build(self):
             doc="Molecular weight kg/mol",
         )
 
+        # Density of water at 25 C
+        self.dens_mass_solvent = Var(
+            within=Reals,
+            initialize=997,
+            units=pyunits.kg / pyunits.m**3,
+            doc="Mass density of water",
+        )
+
         # mass density parameters, eq 4 in Bartholomew
         dens_mass_param_dict = {"0": 995, "1": 756}
         self.dens_mass_param = Var(
@@ -226,7 +234,7 @@ def fix_initialization_states(self):
         # Constraint on water concentration at outlet - unfix in these cases
         for b in self.values():
             if b.config.defined_state is False:
-                b.conc_mol_comp["H2O"].unfix()
+                b.flow_mass_phase_comp["Liq", "H2O"].unfix()
 
     def initialize(
         self,
@@ -491,7 +499,7 @@ def build(self):
         )
 
     # -----------------------------------------------------------------------------
-    # Property Methods
+    # On-demand Property Methods
     def _mass_frac_phase_comp(self):
         self.mass_frac_phase_comp = Var(
             self.params.phase_list,
@@ -788,11 +796,16 @@ def get_enthalpy_flow_terms(self, p):
         return self.enth_flow
 
     # TODO: make property package compatible with dynamics
-    # def get_material_density_terms(self, p, j):
-    #     """Create material density terms."""
-
-    # def get_enthalpy_density_terms(self, p):
-    #     """Create enthalpy density terms."""
+    def get_material_density_terms(self, p, j):
+        """Create material density terms."""
+        if j == "H2O":
+            return self.params.dens_mass_solvent
+        else:
+            return self.conc_mass_phase_comp[p, j]
+
+    def get_energy_density_terms(self, p):
+        """Create enthalpy density terms."""
+        return self.enth_mass_phase[p] * self.dens_mass_phase[p]
 
     def default_material_balance_type(self):
         return MaterialBalanceType.componentTotal

watertap/unit_models/reverse_osmosis_0D.py

@@ -51,8 +51,8 @@ class ReverseOsmosisData(ReverseOsmosisBaseData):
     def _add_feed_side_membrane_channel_and_geometry(self):
         # Build membrane channel control volume
         self.feed_side = MembraneChannel0DBlock(
-            dynamic=False,
-            has_holdup=False,
+            dynamic=self.config.dynamic,
+            has_holdup=self.config.has_holdup,
             property_package=self.config.property_package,
             property_package_args=self.config.property_package_args,
         )

watertap/unit_models/tests/test_cstr.py

@@ -52,6 +52,7 @@
     SingleControlVolumeUnitInitializer,
     InitializationStatus,
 )
+import idaes.logger as idaeslog
 
 # -----------------------------------------------------------------------------
 # Get default solver for testing
@@ -271,13 +272,17 @@ def model(self):
         m.fs.unit.heat_duty[0].fix(0)
         m.fs.unit.deltaP[0].fix(0)
 
+        iscale.calculate_scaling_factors(m)
         return m
 
     @pytest.mark.requires_idaes_solver
     @pytest.mark.component
     def test_general_hierarchical(self, model):
-        initializer = SingleControlVolumeUnitInitializer()
-        initializer.initialize(model.fs.unit)
+        initializer = SingleControlVolumeUnitInitializer(
+            writer_config={"linear_presolve": False}
+        )
+
+        initializer.initialize(model.fs.unit, output_level=idaeslog.DEBUG)
 
         assert initializer.summary[model.fs.unit]["status"] == InitializationStatus.Ok
 

watertap/unit_models/tests/test_reverse_osmosis_0D.py

@@ -9,8 +9,15 @@
 # information, respectively. These files are also available online at the URL
 # "https://github.com/watertap-org/watertap/"
 #################################################################################
-from pyomo.environ import ConcreteModel
+from pyomo.environ import (
+    ConcreteModel,
+    units as pyunits,
+    TransformationFactory,
+    assert_optimal_termination,
+)
 from pyomo.network import Port
+from idaes.core.solvers import petsc
+
 from idaes.core import (
     FlowsheetBlock,
     MaterialBalanceType,
@@ -43,6 +50,10 @@
 from watertap.unit_models.tests.unit_test_harness import UnitTestHarness
 import pytest
 
+from idaes.core.solvers import petsc
+import numpy as np
+
+
 # -----------------------------------------------------------------------------
 # Get default solver for testing
 solver = get_solver()
@@ -115,9 +126,9 @@ def test_default_config_and_build():
         == "fs._time*fs.unit.feed_side.length_domain"
     )
     # test statistics
-    assert number_variables(m) == 134
+    assert number_variables(m) == 135
     assert number_total_constraints(m) == 110
-    assert number_unused_variables(m) == 1
+    assert number_unused_variables(m) == 2
 
 
 def build():
@@ -982,3 +993,104 @@ def configure(self):
         }
 
         return m
+
+
+@pytest.mark.requires_idaes_solver
+@pytest.mark.unit
+def test_RO_dynamic_instantiation():
+    # TODO: add test to check exception for simplest RO0D with dynamics
+
+    m = ConcreteModel()
+    m.fs = FlowsheetBlock(
+        dynamic=True,
+        time_set=list(np.linspace(0, 200, 6)),
+        time_units=pyunits.s,
+    )
+
+    m.fs.properties = props.NaClParameterBlock()
+
+    m.fs.unit = ReverseOsmosis0D(
+        dynamic=True,
+        has_holdup=True,
+        property_package=m.fs.properties,
+        has_pressure_change=True,
+        concentration_polarization_type=ConcentrationPolarizationType.calculated,
+        mass_transfer_coefficient=MassTransferCoefficient.calculated,
+        pressure_change_type=PressureChangeType.calculated,
+        module_type=ModuleType.spiral_wound,
+    )
+
+    time_nfe = len(m.fs.time) - 1
+    TransformationFactory("dae.finite_difference").apply_to(
+        m.fs, nfe=time_nfe, wrt=m.fs.time, scheme="BACKWARD"
+    )
+
+    NaCl_g_per_L_basis = 39
+    NaCl_kg_per_L_basis = NaCl_g_per_L_basis * 1e-3
+    h2o_kg_per_L_basis = 1 - NaCl_kg_per_L_basis  # 0.96
+    desired_feed_flow_start = 22 / 60  # 22 L/min to kg/s
+    NaCl_kg_per_L_start = desired_feed_flow_start / (
+        h2o_kg_per_L_basis / NaCl_kg_per_L_basis + 1
+    )
+    h2o_kg_per_L_start = NaCl_kg_per_L_start * h2o_kg_per_L_basis / NaCl_kg_per_L_basis
+    desired_feed_flow_end = 24 / 60  # 22 L/min to kg/s
+    NaCl_kg_per_L_end = desired_feed_flow_end / (
+        h2o_kg_per_L_basis / NaCl_kg_per_L_basis + 1
+    )
+    h2o_kg_per_L_end = NaCl_kg_per_L_end * h2o_kg_per_L_basis / NaCl_kg_per_L_basis
+    ramp_gradient_NaCl = NaCl_kg_per_L_end - NaCl_kg_per_L_start
+    ramp_gradient_h2o = h2o_kg_per_L_end - h2o_kg_per_L_start
+
+    m.fs.unit.inlet.flow_mass_phase_comp[:, "Liq", "NaCl"].fix(NaCl_kg_per_L_end)
+    m.fs.unit.inlet.flow_mass_phase_comp[:, "Liq", "H2O"].fix(h2o_kg_per_L_end)
+    m.fs.unit.inlet.pressure[:].fix(900 * 6895)  # feed pressure (Pa)
+    ramp_len = 80
+    for i in list([0, 40, ramp_len]):
+        m.fs.unit.inlet.flow_mass_phase_comp[i, "Liq", "NaCl"].fix(NaCl_kg_per_L_start)
+        m.fs.unit.inlet.flow_mass_phase_comp[i, "Liq", "H2O"].fix(h2o_kg_per_L_start)
+        m.fs.unit.inlet.flow_mass_phase_comp[ramp_len + i, "Liq", "NaCl"].fix(
+            NaCl_kg_per_L_start + ramp_gradient_NaCl / ramp_len * i
+        )
+        m.fs.unit.inlet.flow_mass_phase_comp[ramp_len + i, "Liq", "H2O"].fix(
+            h2o_kg_per_L_start + ramp_gradient_h2o / ramp_len * i
+        )
+        m.fs.unit.inlet.pressure[i].fix(800 * 6895)  # feed pressure (Pa)
+        m.fs.unit.inlet.pressure[ramp_len + i].fix(
+            (800 + 100 / ramp_len * i) * 6895
+        )  # feed pressure (Pa)
+
+    m.fs.unit.inlet.temperature[:].fix(293.15)  # feed temperature (K)
+
+    m.fs.unit.area.fix(7.4)  # membrane area (m^2)
+    m.fs.unit.A_comp.fix(3.75e-12)  # membrane water permeability (m/Pa/s)
+    m.fs.unit.B_comp.fix(3e-8)  # membrane salt permeability (m/s)
+    m.fs.unit.permeate.pressure[:].fix(101325)  # permeate pressure (Pa)
+
+    m.fs.unit.feed_side.channel_height.fix(0.001)
+    m.fs.unit.feed_side.spacer_porosity.fix(0.729)  # 72.9%
+    m.fs.unit.length.fix(1.016 - 2 * 0.0267)  # m
+
+    m.fs.unit.feed_side.material_accumulation[:, :, :].value = 0
+    m.fs.unit.feed_side.material_accumulation[0, :, :].fix(0)
+
+    assert not hasattr(m.fs.unit.feed_side, "energy_accumulation")
+
+    m.fs.properties.set_default_scaling("flow_mass_phase_comp", 1, index=("Liq", "H2O"))
+    m.fs.properties.set_default_scaling(
+        "flow_mass_phase_comp", 1e2, index=("Liq", "NaCl")
+    )
+
+    iscale.set_scaling_factor(m.fs.unit.area, 1e-2)
+
+    iscale.calculate_scaling_factors(m)
+
+    m.fs.unit.initialize()
+
+    iscale.calculate_scaling_factors(m)
+
+    results = petsc.petsc_dae_by_time_element(
+        m,
+        time=m.fs.time,
+    )
+    for result in results.results:
+        assert_optimal_termination(result)