Enable RO0D to be dynamic (#1471)

* revised dyanmic and has holdup for membrane channel base/RO 0D, but prop pack needs get_material_density_terms * Changed dynamic and has_holdup to not be hardcoded to False * run black * cleanup test * Fixed NaCl property package as per TIm's instructions to allow dynamics * add material_density by component * blk * put mat dens expr in right place * fix mat dens epr brackets * fix ro test * blk * fix unintentional addition of conc_mass_phase_comp to state vars * remove commented out code * add new instance with simple configs * black * Added numbers to initialize RO0D unit_models * raise config error for simple RO * Added pressure change back in to be identical to Adam's specs * setup.py point to dynamic_cv0d * add back pyomo * black * revert debug changes after idaes PR fix; add check to tests to check units on material_accumulation * cleanup and check material_holdup * adjust test * blk * initialize without the assert_units_consistent check in place * mixed_permeate as CV0D * black RO-base * add permeate port to mixed_permeate * revert debug changes after idaes PR fix; add check to tests to check units on material_accumulation * cleanup and check material_holdup * adjust test * blk * initialize without the assert_units_consistent check in place * undo RO0D and RObase changes * revert setup.py * black again * add idaes temporary dependency * black * add idaes dependency with CV0D fix * blk * move assertion * Currently DOF=2 * petsc succeeds with DOF=0 * cleanup * bring my changes back in after merging with fahim's last commits * remove unused imports * one more import to delete * add petsc to github actions * skip test for workflows that don't have petsc available; note, tests work and pass for windows/linux checks * Nonzero accumulation. Can get dynamic change now. * Getting a plot for RO0D now. * Plotting permeate NaCl conc. * Plotting feed outlet C_NaCl * Blacked RO0D test * Plots 5 things. * Preliminary results obtained. * try to fix failed mac test and run black * revised dyanmic and has holdup for membrane channel base/RO 0D, but prop pack needs get_material_density_terms * add material_density by component * fix mat dens epr brackets * add new instance with simple configs * raise config error for simple RO * black * revert debug changes after idaes PR fix; add check to tests to check units on material_accumulation * revert debug changes after idaes PR fix; add check to tests to check units on material_accumulation * petsc succeeds with DOF=0 * Currently DOF=2 * bring my changes back in after merging with fahim's last commits * remove unused imports * Validate dynamic results. * Validate dynamic model. * Validate dynamic model. * black * try to fix test_cstr * bk * test commit * black again * remove unused import * Revert "try to fix failed mac test and run black" This reverts commit 546d4d7. * test ro back up and running * run black and remove extra import * cleanup RO base * Apply suggestions from code review Co-authored-by: Ludovico Bianchi <[email protected]> * Update setup.py Co-authored-by: Ludovico Bianchi <[email protected]> * fix checks * probe failure * bring back reverted commit * Revert "Merge branch 'main' into dynamic-exploratory2" This reverts commit 1cdd1a2, reversing changes made to a282f18. * remove file that surprisingly made its way into a commit * point to 2.6rc * adjust HRT constraint on cstr and check failing macos test * remove unused import * black * reset membrance_channel0d to match main * merge main flowsheets * bring back other changes that are in main already * revert hydraulic retention time constraint * checkpoint: add scaling to cstr but didn't note any affect on BSM2-P solver iterations; also replacing aeration tank on bsm2-p instead of cstr_injection with flowsheet level vars and constraints * add requires idaes solver mark for dynamic ro test * cleanup * Revert "checkpoint: add scaling to cstr but didn't note any affect on BSM2-P solver iterations; also replacing aeration tank on bsm2-p instead of cstr_injection with flowsheet level vars and constraints" This reverts commit bd2caa8. * black * address pylint --------- Co-authored-by: Adam Atia <[email protected]> Co-authored-by: kurbansitterley <[email protected]> Co-authored-by: Ludovico Bianchi <[email protected]>
watertap-org · Oct 30, 2024 · 53bd64e · 53bd64e
1 parent 8c651de
commit 53bd64e
Show file tree

Hide file tree

Showing 6 changed files with 160 additions and 20 deletions.
diff --git a/.github/workflows/checks.yml b/.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: |

diff --git a/watertap/core/membrane_channel_base.py b/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

diff --git a/watertap/property_models/NaCl_prop_pack.py b/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

diff --git a/watertap/unit_models/reverse_osmosis_0D.py b/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,
         )

diff --git a/watertap/unit_models/tests/test_cstr.py b/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
 

diff --git a/watertap/unit_models/tests/test_reverse_osmosis_0D.py b/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)