diff --git a/.github/workflows/checks.yml b/.github/workflows/checks.yml index cfb73390f3..585d080153 100644 --- 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 index 7f0aeb050f..0b3f365599 100644 --- 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 index 56deaf5849..31f0b0157e 100644 --- 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 index 235be5d762..e5b8342075 100644 --- 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 index 41c68ef585..ba53472386 100644 --- 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 index daa994dcfe..ee10586c2a 100644 --- 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)