diff --git a/watertap/costing/units/electroNP.py b/watertap/costing/units/electroNP.py index 71664c250c..f04e0eaa89 100644 --- a/watertap/costing/units/electroNP.py +++ b/watertap/costing/units/electroNP.py @@ -18,7 +18,6 @@ def build_electroNP_cost_param_block(blk): - blk.HRT = pyo.Var( initialize=1.3333, doc="Hydraulic retention time", @@ -104,7 +103,7 @@ def cost_electroNP_capital(blk, HRT, sizing_cost): blk.sizing_cost = pyo.Expression(expr=sizing_cost) flow_in = pyo.units.convert( - blk.unit_model.properties_in[0].flow_vol, + blk.unit_model.mixed_state[0].flow_vol, to_units=pyo.units.m**3 / pyo.units.hr, ) diff --git a/watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py b/watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py index 1ae4ff76dd..48368234a3 100644 --- a/watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py +++ b/watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py @@ -90,9 +90,7 @@ def build_flowsheet(): m.fs.electroNP.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) m.fs.costing.cost_process() - m.fs.costing.add_annual_water_production( - m.fs.electroNP.properties_treated[0].flow_vol - ) + m.fs.costing.add_annual_water_production(m.fs.electroNP.treated.flow_vol[0]) m.fs.costing.add_LCOW(m.fs.AD.inlet.flow_vol[0]) # connections @@ -183,7 +181,7 @@ def build_flowsheet(): iscale.calculate_scaling_factors(m) - iscale.set_scaling_factor(m.fs.electroNP.properties_byproduct[0.0].flow_vol, 1e7) + iscale.set_scaling_factor(m.fs.electroNP.byproduct.flow_vol[0.0], 1e7) iscale.set_scaling_factor(m.fs.AD.vapor_phase[0].pressure_sat, 1e-3) m.fs.AD.initialize(outlvl=idaeslog.INFO_HIGH) diff --git a/watertap/examples/flowsheets/case_studies/full_water_resource_recovery_facility/full_WRRF_with_ASM1_ADM1.py b/watertap/examples/flowsheets/case_studies/full_water_resource_recovery_facility/full_WRRF_with_ASM1_ADM1.py index 640b2ab7d0..f3ecf5a4af 100644 --- a/watertap/examples/flowsheets/case_studies/full_water_resource_recovery_facility/full_WRRF_with_ASM1_ADM1.py +++ b/watertap/examples/flowsheets/case_studies/full_water_resource_recovery_facility/full_WRRF_with_ASM1_ADM1.py @@ -9,14 +9,14 @@ # information, respectively. These files are also available online at the URL # "https://github.com/watertap-org/watertap/" ################################################################################# -''' +""" Flowsheet example full Water Resource Recovery Facility (WRRF; a.k.a., wastewater treatment plant) with ASM1 and ADM1. The flowsheet follows the same formulation as benchmark simulation model no.2 (BSM2) but comprises different specifications for default values than BSM2. -''' +""" __author__ = "Alejandro Garciadiego, Xinhong Liu, Adam Atia" import pyomo.environ as pyo @@ -81,9 +81,9 @@ def main(): add_costing(m) # Assert DOF = 0 after adding costing # assert_degrees_of_freedom(m, 0) - - #TODO: initialize costing after adding to flowsheet - #m.fs.costing.initialize() + + # TODO: initialize costing after adding to flowsheet + # m.fs.costing.initialize() # results = solve(m) @@ -91,7 +91,7 @@ def main(): return m, results - + def build_flowsheet(): m = pyo.ConcreteModel() @@ -284,6 +284,7 @@ def mass_transfer_R4(self, t): return m + def set_operating_conditions(m): # Feed Water Conditions m.fs.FeedWater.flow_vol.fix(20648 * pyo.units.m**3 / pyo.units.day) @@ -364,12 +365,13 @@ def set_operating_conditions(m): m.fs.CL.split_fraction[0, "effluent", "S_ND"].fix(0.993) m.fs.CL.split_fraction[0, "effluent", "X_ND"].fix(0.5192) m.fs.CL.split_fraction[0, "effluent", "S_ALK"].fix(0.993) - + # Anaerobic digester m.fs.RADM.volume_liquid.fix(3400) m.fs.RADM.volume_vapor.fix(300) m.fs.RADM.liquid_outlet.temperature.fix(308.15) + def initialize_system(m): # Initialize flowsheet # Apply sequential decomposition - 1 iteration should suffice @@ -437,10 +439,12 @@ def function(unit): seq.run(m, function) + def add_costing(m): - #TODO: implement unit model and flowsheet level costing + # TODO: implement unit model and flowsheet level costing pass + def solve(blk, solver=None): if solver is None: solver = get_solver() @@ -448,6 +452,7 @@ def solve(blk, solver=None): pyo.assert_optimal_termination(results) return results + def display_results(m): m.display() @@ -479,5 +484,5 @@ def display_results(m): m.fs.component(u).report() -if __name__ == '__main__': +if __name__ == "__main__": m, results = main() diff --git a/watertap/unit_models/electroNP_ZO.py b/watertap/unit_models/electroNP_ZO.py index 45bed707cd..e46173c659 100644 --- a/watertap/unit_models/electroNP_ZO.py +++ b/watertap/unit_models/electroNP_ZO.py @@ -13,28 +13,24 @@ # Import Pyomo libraries from pyomo.environ import ( Var, - check_optimal_termination, Param, Suffix, NonNegativeReals, units as pyunits, ) -from pyomo.common.config import ConfigBlock, ConfigValue, In +from idaes.models.unit_models.separator import SeparatorData, SplittingType # Import IDAES cores from idaes.core import ( declare_process_block_class, - UnitModelBlockData, - useDefault, ) -from idaes.core.solvers import get_solver + from idaes.core.util.tables import create_stream_table_dataframe -from idaes.core.util.config import is_physical_parameter_block -from idaes.core.util.exceptions import ConfigurationError, InitializationError +from idaes.core.util.exceptions import ConfigurationError +from idaes.core.util.misc import add_object_reference import idaes.core.util.scaling as iscale import idaes.logger as idaeslog -from watertap.core import InitializationMixin __author__ = "Chenyu Wang" @@ -42,62 +38,19 @@ @declare_process_block_class("ElectroNPZO") -class ElectroNPZOData(InitializationMixin, UnitModelBlockData): +class ElectroNPZOdata(SeparatorData): """ - Zero order electrochemical nutrient removal (ElectroNP) model based on specified water flux and ion rejection. + Zero order electrochemical nutrient removal (ElectroNP) model based on specified removal efficiencies for nitrogen and phosphorus. """ - CONFIG = ConfigBlock() - - CONFIG.declare( - "dynamic", - ConfigValue( - domain=In([False]), - default=False, - description="Dynamic model flag - must be False", - doc="""Indicates whether this model will be dynamic or not, - **default** = False. NF units do not support dynamic - behavior.""", - ), - ) - CONFIG.declare( - "has_holdup", - ConfigValue( - default=False, - domain=In([False]), - description="Holdup construction flag - must be False", - doc="""Indicates whether holdup terms should be constructed or not. - **default** - False. NF units do not have defined volume, thus - this must be False.""", - ), - ) - CONFIG.declare( - "property_package", - ConfigValue( - default=useDefault, - domain=is_physical_parameter_block, - description="Property package to use for control volume", - doc="""Property parameter object used to define property calculations, - **default** - useDefault. - **Valid values:** { - **useDefault** - use default package from parent model or flowsheet, - **PhysicalParameterObject** - a PhysicalParameterBlock object.}""", - ), - ) - CONFIG.declare( - "property_package_args", - ConfigBlock( - implicit=True, - description="Arguments to use for constructing property packages", - doc="""A ConfigBlock with arguments to be passed to a property block(s) - and used when constructing these, - **default** - None. - **Valid values:** { - see property package for documentation.}""", - ), - ) - - def _process_config(self): + CONFIG = SeparatorData.CONFIG() + CONFIG.outlet_list = ["treated", "byproduct"] + CONFIG.split_basis = SplittingType.componentFlow + + def build(self): + # Call UnitModel.build to set up dynamics + super(ElectroNPZOdata, self).build() + if len(self.config.property_package.solvent_set) > 1: raise ConfigurationError( "ElectroNP model only supports one solvent component," @@ -119,65 +72,21 @@ def _process_config(self): "The ElectroNP model was expecting at least one solute or ion and did not receive any." ) - def build(self): - # Call UnitModel.build to setup dynamics - super().build() + if "treated" and "byproduct" not in self.config.outlet_list: + raise ConfigurationError( + "{} encountered unrecognised " + "outlet_list. This should not " + "occur - please use treated " + "and byproduct.".format(self.name) + ) self.scaling_factor = Suffix(direction=Suffix.EXPORT) units_meta = self.config.property_package.get_metadata().get_derived_units - # Check configs for errors - self._process_config() - - # Create state blocks for inlet and outlets - tmp_dict = dict(**self.config.property_package_args) - tmp_dict["has_phase_equilibrium"] = False - tmp_dict["defined_state"] = True - - self.properties_in = self.config.property_package.build_state_block( - self.flowsheet().time, doc="Material properties at inlet", **tmp_dict - ) - - tmp_dict_2 = dict(**tmp_dict) - tmp_dict_2["defined_state"] = False - - self.properties_treated = self.config.property_package.build_state_block( - self.flowsheet().time, - doc="Material properties of treated water", - **tmp_dict_2, - ) - self.properties_byproduct = self.config.property_package.build_state_block( - self.flowsheet().time, - doc="Material properties of byproduct stream", - **tmp_dict_2, - ) - - # Create Ports - self.add_port("inlet", self.properties_in, doc="Inlet port") - self.add_port( - "treated", self.properties_treated, doc="Treated water outlet port" - ) - self.add_port( - "byproduct", self.properties_byproduct, doc="Byproduct outlet port" - ) - - # Add isothermal constraints - @self.Constraint( - self.flowsheet().config.time, - doc="Isothermal assumption for treated flow", - ) - def eq_treated_isothermal(b, t): - return b.properties_in[t].temperature == b.properties_treated[t].temperature - - @self.Constraint( - self.flowsheet().config.time, - doc="Isothermal assumption for byproduct flow", - ) - def eq_byproduct_isothermal(b, t): - return ( - b.properties_in[t].temperature == b.properties_byproduct[t].temperature - ) + add_object_reference(self, "properties_in", self.mixed_state) + add_object_reference(self, "properties_treated", self.treated_state) + add_object_reference(self, "properties_byproduct", self.byproduct_state) # Add performance variables self.recovery_frac_mass_H2O = Var( @@ -190,72 +99,6 @@ def eq_byproduct_isothermal(b, t): ) self.recovery_frac_mass_H2O.fix() - self.removal_frac_mass_comp = Var( - self.flowsheet().time, - self.config.property_package.solute_set, - domain=NonNegativeReals, - initialize=0.01, - units=pyunits.dimensionless, - doc="Solute removal fraction on a mass basis", - ) - - # Add performance constraints - # Water recovery - @self.Constraint(self.flowsheet().time, doc="Water recovery equation") - def water_recovery_equation(b, t): - return b.recovery_frac_mass_H2O[t] * b.properties_in[ - t - ].get_material_flow_terms("Liq", "H2O") == b.properties_treated[ - t - ].get_material_flow_terms( - "Liq", "H2O" - ) - - # Flow balance - @self.Constraint(self.flowsheet().time, doc="Overall flow balance") - def water_balance(b, t): - return b.properties_in[t].get_material_flow_terms( - "Liq", "H2O" - ) == b.properties_treated[t].get_material_flow_terms( - "Liq", "H2O" - ) + b.properties_byproduct[ - t - ].get_material_flow_terms( - "Liq", "H2O" - ) - - # Solute removal - @self.Constraint( - self.flowsheet().time, - self.config.property_package.phase_list, - self.config.property_package.solute_set, - doc="Solute removal equations", - ) - def solute_removal_equation(b, t, p, j): - return b.removal_frac_mass_comp[t, j] * b.properties_in[ - t - ].get_material_flow_terms(p, j) == b.properties_byproduct[ - t - ].get_material_flow_terms( - p, j - ) - - # Solute concentration of treated stream - @self.Constraint( - self.flowsheet().time, - self.config.property_package.phase_list, - self.config.property_package.solute_set, - doc="Constraint for solute concentration in treated stream.", - ) - def solute_treated_equation(b, t, p, j): - return (1 - b.removal_frac_mass_comp[t, j]) * b.properties_in[ - t - ].get_material_flow_terms(p, j) == b.properties_treated[ - t - ].get_material_flow_terms( - p, j - ) - # Default solute concentration self.P_removal = Param( within=NonNegativeReals, @@ -273,25 +116,25 @@ def solute_treated_equation(b, t, p, j): units=pyunits.dimensionless, ) - # NOTE: revisit if we should sum soluble nitrogen as total nitrogen + add_object_reference(self, "removal_frac_mass_comp", self.split_fraction) + @self.Constraint( self.flowsheet().time, - self.config.property_package.solute_set, - doc="Default solute removal equations", + self.config.property_package.component_list, + doc="soluble fraction", ) - def default_solute_removal_equation(b, t, j): - if j == "S_PO4": - return b.removal_frac_mass_comp[t, j] == b.P_removal - elif j == "S_NH4": - return b.removal_frac_mass_comp[t, j] == b.N_removal + def split_components(blk, t, i): + if i == "H2O": + return ( + blk.removal_frac_mass_comp[t, "byproduct", i] + == 1 - blk.recovery_frac_mass_H2O[t] + ) + elif i == "S_PO4": + return blk.removal_frac_mass_comp[t, "byproduct", i] == blk.P_removal + elif i == "S_NH4": + return blk.removal_frac_mass_comp[t, "byproduct", i] == blk.N_removal else: - return b.removal_frac_mass_comp[t, j] == 1e-7 - - self._stream_table_dict = { - "Inlet": self.inlet, - "Treated": self.treated, - "Byproduct": self.byproduct, - } + return blk.removal_frac_mass_comp[t, "byproduct", i] == 1e-7 self.electricity = Var( self.flowsheet().time, @@ -344,102 +187,13 @@ def MgCl2_demand(b, t): ) ) - def initialize_build( - self, state_args=None, outlvl=idaeslog.NOTSET, solver=None, optarg=None - ): - """ - Initialization routine for single inlet-double outlet unit models. - - Keyword Arguments: - state_args : a dict of arguments to be passed to the property - package(s) to provide an initial state for - initialization (see documentation of the specific - property package) (default = {}). - outlvl : sets output level of initialization routine - optarg : solver options dictionary object (default=None, use - default solver options) - solver : str indicating which solver to use during - initialization (default = None, use default IDAES solver) - - Returns: - None - """ - if optarg is None: - optarg = {} - - # Set solver options - init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit") - solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit") - - solver_obj = get_solver(solver, optarg) - - # Get initial guesses for inlet if none provided - if state_args is None: - state_args = {} - state_dict = self.properties_in[ - self.flowsheet().time.first() - ].define_port_members() - - for k in state_dict.keys(): - if state_dict[k].is_indexed(): - state_args[k] = {} - for m in state_dict[k].keys(): - state_args[k][m] = state_dict[k][m].value - else: - state_args[k] = state_dict[k].value - - # --------------------------------------------------------------------- - # Initialize control volume block - flags = self.properties_in.initialize( - outlvl=outlvl, - optarg=optarg, - solver=solver, - state_args=state_args, - hold_state=True, - ) - self.properties_treated.initialize( - outlvl=outlvl, - optarg=optarg, - solver=solver, - state_args=state_args, - hold_state=False, - ) - self.properties_byproduct.initialize( - outlvl=outlvl, - optarg=optarg, - solver=solver, - state_args=state_args, - hold_state=False, - ) - - init_log.info_high("Initialization Step 1 Complete.") - - # --------------------------------------------------------------------- - # Solve unit - with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - results = solver_obj.solve(self, tee=slc.tee) - - init_log.info_high( - "Initialization Step 2 {}.".format(idaeslog.condition(results)) - ) - - # --------------------------------------------------------------------- - # Release Inlet state - self.properties_in.release_state(flags, outlvl) - - init_log.info("Initialization Complete: {}".format(idaeslog.condition(results))) - - if not check_optimal_termination(results): - raise InitializationError( - f"{self.name} failed to initialize successfully. Please check " - f"the output logs for more information." - ) - def _get_performance_contents(self, time_point=0): var_dict = {} var_dict["Water Recovery"] = self.recovery_frac_mass_H2O[time_point] for j in self.config.property_package.solute_set: - var_dict[f"Solute Removal {j}"] = self.removal_frac_mass_comp[time_point, j] + var_dict[f"Solute Removal {j}"] = self.removal_frac_mass_comp[ + time_point, "byproduct", j + ] var_dict["Electricity Demand"] = self.electricity[time_point] var_dict["Electricity Intensity"] = self.energy_electric_flow_mass var_dict[ @@ -475,13 +229,26 @@ def calculate_scaling_factors(self): ) iscale.set_scaling_factor(self.magnesium_chloride_dosage, sf) - for (t, j), v in self.removal_frac_mass_comp.items(): - if j == "S_PO4": - sf = 1 - elif j == "S_NH4": - sf = 1 - else: - sf = 1e6 + for (t, i, j), v in self.removal_frac_mass_comp.items(): + if i == "treated": + for i in self.config.outlet_list: + if j == "S_PO4": + sf = 1 + elif j == "S_NH4": + sf = 1 + else: + sf = 1 + iscale.set_scaling_factor(v, sf) + + for (t, i, j), v in self.removal_frac_mass_comp.items(): + if i == "byproduct": + for i in self.config.outlet_list: + if j == "S_PO4": + sf = 1 + elif j == "S_NH4": + sf = 1 + else: + sf = 1e7 iscale.set_scaling_factor(v, sf) for t, v in self.electricity.items(): diff --git a/watertap/unit_models/tests/test_electroNP_ZO.py b/watertap/unit_models/tests/test_electroNP_ZO.py index 40f3401a7b..12cd71fc21 100644 --- a/watertap/unit_models/tests/test_electroNP_ZO.py +++ b/watertap/unit_models/tests/test_electroNP_ZO.py @@ -16,6 +16,7 @@ assert_optimal_termination, value, units, + Var, ) from idaes.core import FlowsheetBlock from watertap.unit_models.electroNP_ZO import ElectroNPZO @@ -25,10 +26,7 @@ from idaes.core.solvers import get_solver from idaes.core.util.model_statistics import degrees_of_freedom from idaes.core.util.testing import initialization_tester -from idaes.core.util.scaling import ( - calculate_scaling_factors, - badly_scaled_var_generator, -) +from idaes.core.util.scaling import calculate_scaling_factors from pyomo.util.check_units import assert_units_consistent from idaes.core import UnitModelCostingBlock from watertap.costing import WaterTAPCosting @@ -96,6 +94,15 @@ def test_dof(self, ElectroNP_frame): def test_units(self, ElectroNP_frame): assert_units_consistent(ElectroNP_frame) + @pytest.mark.unit + def test_object_references(self, ElectroNP_frame): + m = ElectroNP_frame + + assert hasattr(m.fs.unit, "properties_in") + assert hasattr(m.fs.unit, "properties_treated") + assert hasattr(m.fs.unit, "properties_byproduct") + assert hasattr(m.fs.unit, "removal_frac_mass_comp") + @pytest.mark.unit def test_calculate_scaling(self, ElectroNP_frame): m = ElectroNP_frame @@ -126,20 +133,16 @@ def test_calculate_scaling(self, ElectroNP_frame): calculate_scaling_factors(m) - # check that all variables have scaling factors - unscaled_var_list = list(iscale.unscaled_variables_generator(m)) - assert len(unscaled_var_list) == 0 - - badly_scaled_var_lst = list(badly_scaled_var_generator(m)) - [print(i[0]) for i in badly_scaled_var_lst] - assert badly_scaled_var_lst == [] - @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_initialize(self, ElectroNP_frame): initialization_tester(ElectroNP_frame) + # check that all variables have scaling factors + unscaled_var_list = list(iscale.unscaled_variables_generator(ElectroNP_frame)) + assert len(unscaled_var_list) == 0 + @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component