diff --git a/idaes/models/properties/modular_properties/base/generic_property.py b/idaes/models/properties/modular_properties/base/generic_property.py index bef5ddd22c..726dd702dd 100644 --- a/idaes/models/properties/modular_properties/base/generic_property.py +++ b/idaes/models/properties/modular_properties/base/generic_property.py @@ -2726,7 +2726,10 @@ def initialize( # For systems where the state variables fully define the # phase equilibrium, we cannot activate the equilibrium # constraint at this stage. - if "flow_mol_phase_comp" not in b.define_state_vars(): + if ( + "flow_mol_phase_comp" not in b.define_state_vars() + and "mole_frac_phase_comp" not in b.define_state_vars() + ): c.activate() for pp in b.params._pe_pairs: @@ -2858,21 +2861,27 @@ def build(self): self.params.config.state_definition.define_state(self) # Add equilibrium temperature variable if required - if self.params.config.phases_in_equilibrium is not None and ( - not self.config.defined_state or self.always_flash - ): - - t_units = self.params.get_metadata().default_units.TEMPERATURE - if self.temperature.value is not None: - t_value = value(self.temperature) + if self.params.config.phases_in_equilibrium is not None: + if self.always_flash: + has_phase_equilibrium = True + elif self.config.defined_state: + has_phase_equilibrium = False else: - t_value = None - self._teq = Var( - self.params._pe_pairs, - initialize=t_value, - doc="Temperature for calculating phase equilibrium", - units=t_units, - ) + has_phase_equilibrium = self.config.has_phase_equilibrium + + if has_phase_equilibrium: + t_units = self.params.get_metadata().default_units.TEMPERATURE + if self.temperature.value is not None: + t_value = value(self.temperature) + else: + t_value = None + self._teq = Var( + self.params._pe_pairs, + initialize=t_value, + doc="Temperature for calculating phase equilibrium", + units=t_units, + bounds=self.temperature.bounds, + ) # Create common components for each property package for p in self.phase_list: @@ -2888,37 +2897,42 @@ def enth_mol_eqn(b): b.enth_mol_phase[p] * b.phase_frac[p] for p in b.phase_list ) - # Add phase equilibrium constraints if necessary - if self.params.config.phases_in_equilibrium is not None and ( - not self.config.defined_state or self.always_flash - ): + if self.params.config.phases_in_equilibrium is not None: + if self.always_flash: + has_phase_equilibrium = True + elif self.config.defined_state: + has_phase_equilibrium = False + else: + has_phase_equilibrium = self.config.has_phase_equilibrium - pe_form_config = self.params.config.phase_equilibrium_state - for pp in self.params._pe_pairs: - pe_form_config[pp].phase_equil(self, pp) - - def rule_equilibrium(b, phase1, phase2, j): - if (phase1, j) not in b.phase_component_set or ( - phase2, - j, - ) not in b.phase_component_set: - return Constraint.Skip - config = b.params.get_component(j).config - try: - e_mthd = config.phase_equilibrium_form[ - (phase1, phase2) - ].return_expression - except KeyError: - e_mthd = config.phase_equilibrium_form[ - (phase2, phase1) - ].return_expression - if e_mthd is None: - raise GenericPropertyPackageError(b, "phase_equilibrium_form") - return e_mthd(self, phase1, phase2, j) + if has_phase_equilibrium: - self.equilibrium_constraint = Constraint( - self.params._pe_pairs, self.component_list, rule=rule_equilibrium - ) + pe_form_config = self.params.config.phase_equilibrium_state + for pp in self.params._pe_pairs: + pe_form_config[pp].phase_equil(self, pp) + + def rule_equilibrium(b, phase1, phase2, j): + if (phase1, j) not in b.phase_component_set or ( + phase2, + j, + ) not in b.phase_component_set: + return Constraint.Skip + config = b.params.get_component(j).config + try: + e_mthd = config.phase_equilibrium_form[ + (phase1, phase2) + ].return_expression + except KeyError: + e_mthd = config.phase_equilibrium_form[ + (phase2, phase1) + ].return_expression + if e_mthd is None: + raise GenericPropertyPackageError(b, "phase_equilibrium_form") + return e_mthd(self, phase1, phase2, j) + + self.equilibrium_constraint = Constraint( + self.params._pe_pairs, self.component_list, rule=rule_equilibrium + ) # Add inherent reaction constraints if necessary if self.params.has_inherent_reactions and ( @@ -5503,12 +5517,15 @@ def _temperature_pressure_bubble_dew(b, name): abbrv = "t" + abbrv bounds = (b.temperature.lb, b.temperature.ub) units = b.params.get_metadata().default_units.TEMPERATURE + init = b.temperature.value elif splt[0] == "pressure": abbrv = "p" + abbrv bounds = (b.pressure.lb, b.pressure.ub) units_meta = b.params.get_metadata().derived_units units = units_meta.PRESSURE + init = b.pressure.value else: + init = None _raise_dev_burnt_toast() docstring_var += splt[0] @@ -5521,7 +5538,13 @@ def _temperature_pressure_bubble_dew(b, name): setattr( b, name, - Var(b.params._pe_pairs, doc=docstring_var, bounds=bounds, units=units), + Var( + b.params._pe_pairs, + doc=docstring_var, + bounds=bounds, + initialize=init, + units=units, + ), ) setattr( b, diff --git a/idaes/models/properties/modular_properties/examples/tests/test_BTIdeal_FpTPxpc.py b/idaes/models/properties/modular_properties/examples/tests/test_BTIdeal_FpTPxpc.py new file mode 100644 index 0000000000..9eaef6cc98 --- /dev/null +++ b/idaes/models/properties/modular_properties/examples/tests/test_BTIdeal_FpTPxpc.py @@ -0,0 +1,271 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2024 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Tests for feed with flash with the FpTPxpc state definition. +Authors: Tanner Polley +""" + +import pytest + +from pyomo.environ import ( + check_optimal_termination, + ConcreteModel, + Set, + value, + Var, + units as pyunits, +) +from pyomo.util.check_units import assert_units_consistent +from pyomo.common.unittest import assertStructuredAlmostEqual + +from idaes.core import FlowsheetBlock, Component +from idaes.models.unit_models.feed_flash import FeedFlash +from idaes.models.properties.activity_coeff_models.BTX_activity_coeff_VLE import ( + BTXParameterBlock, +) +from idaes.models.properties.modular_properties import GenericParameterBlock +from idaes.models.properties.modular_properties.state_definitions import FpTPxpc +from idaes.models.properties.modular_properties.examples.BT_ideal import ( + configuration, +) + +from idaes.models.properties.modular_properties.phase_equil import SmoothVLE +from idaes.models.properties.tests.test_harness import PropertyTestHarness + +from idaes.core.util.testing import PhysicalParameterTestBlock, initialization_tester +from idaes.core.solvers import get_solver + +from idaes.core.util import DiagnosticsToolbox + +from copy import deepcopy + + +# ----------------------------------------------------------------------------- +# Get default solver for testing +solver = get_solver("ipopt_v2") + + +# ----------------------------------------------------------------------------- + +configuration = deepcopy(configuration) + +configuration["state_definition"] = FpTPxpc +configuration["state_bounds"] = { + "flow_mol_phase": (0, 100, 1000, pyunits.mol / pyunits.s), + "temperature": (273.15, 300, 450, pyunits.K), + "pressure": (5e4, 1e5, 1e6, pyunits.Pa), +} + + +def _as_quantity(x): + unit = pyunits.get_units(x) + if unit is None: + unit = pyunits.dimensionless + return value(x) * unit._get_pint_unit() + + +class TestBTIdeal(PropertyTestHarness): + def configure(self): + self.prop_pack = GenericParameterBlock + self.param_args = configuration + self.prop_args = {} + self.has_density_terms = True + + +class TestParamBlock(object): + @pytest.mark.unit + def test_build(self): + model = ConcreteModel() + model.params = GenericParameterBlock(**configuration) + + assert isinstance(model.params.phase_list, Set) + assert len(model.params.phase_list) == 2 + for i in model.params.phase_list: + assert i in ["Liq", "Vap"] + assert model.params.Liq.is_liquid_phase() + assert model.params.Vap.is_vapor_phase() + + assert isinstance(model.params.component_list, Set) + assert len(model.params.component_list) == 2 + for i in model.params.component_list: + assert i in ["benzene", "toluene"] + assert isinstance(model.params.get_component(i), Component) + + assert isinstance(model.params._phase_component_set, Set) + assert len(model.params._phase_component_set) == 4 + for i in model.params._phase_component_set: + assert i in [ + ("Liq", "benzene"), + ("Liq", "toluene"), + ("Vap", "benzene"), + ("Vap", "toluene"), + ] + + assert model.params.config.state_definition == FpTPxpc + + assertStructuredAlmostEqual( + model.params.config.state_bounds, + { + "flow_mol_phase": (0, 100, 1000, pyunits.mol / pyunits.s), + "temperature": (273.15, 300, 450, pyunits.K), + "pressure": (5e4, 1e5, 1e6, pyunits.Pa), + }, + item_callback=_as_quantity, + ) + + assert model.params.config.phase_equilibrium_state == { + ("Vap", "Liq"): SmoothVLE + } + + assert isinstance(model.params.phase_equilibrium_idx, Set) + assert len(model.params.phase_equilibrium_idx) == 2 + for i in model.params.phase_equilibrium_idx: + assert i in ["PE1", "PE2"] + + assert model.params.phase_equilibrium_list == { + "PE1": ["benzene", ("Vap", "Liq")], + "PE2": ["toluene", ("Vap", "Liq")], + } + + assert model.params.pressure_ref.value == 1e5 + assert model.params.temperature_ref.value == 300 + + assert model.params.benzene.mw.value == 78.1136e-3 + assert model.params.benzene.pressure_crit.value == 48.9e5 + assert model.params.benzene.temperature_crit.value == 562.2 + + assert model.params.toluene.mw.value == 92.1405e-3 + assert model.params.toluene.pressure_crit.value == 41e5 + assert model.params.toluene.temperature_crit.value == 591.8 + + assert_units_consistent(model) + + +# ----------------------------------------------------------------------------- +class TestBTIdeal_FpTPxpc(object): + @pytest.fixture(scope="class") + def model(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = GenericParameterBlock(**configuration) + + m.fs.unit = FeedFlash(property_package=m.fs.properties) + + m.fs.unit.temperature.fix(368.0) + m.fs.unit.pressure.fix(101325.0) + + m.fs.unit.flow_mol_phase[0, "Liq"].fix(0.5) + m.fs.unit.flow_mol_phase[0, "Vap"].fix(0.5) + m.fs.unit.mole_frac_phase_comp[0, "Liq", "benzene"].fix(0.5) + m.fs.unit.mole_frac_phase_comp[0, "Liq", "toluene"].fix(0.5) + m.fs.unit.mole_frac_phase_comp[0, "Vap", "benzene"].fix(0.5) + m.fs.unit.mole_frac_phase_comp[0, "Vap", "toluene"].fix(0.5) + + return m + + @pytest.mark.build + @pytest.mark.unit + def test_build(self, model): + + # Check state variable values and bounds + assert isinstance(model.fs.unit.flow_mol_phase, Var) + for p in model.fs.unit.flow_mol_phase: + assert value(model.fs.unit.flow_mol_phase[p]) == 0.5 + assert model.fs.unit.flow_mol_phase[p].ub == 1000 + assert model.fs.unit.flow_mol_phase[p].lb == 0 + + assert isinstance(model.fs.unit.pressure, Var) + assert value(model.fs.unit.pressure[0]) == 101325 + assert model.fs.unit.pressure[0].ub == 1e6 + assert model.fs.unit.pressure[0].lb == 5e4 + + assert isinstance(model.fs.unit.temperature, Var) + assert value(model.fs.unit.temperature[0]) == 368 + assert model.fs.unit.temperature[0].ub == 450 + assert model.fs.unit.temperature[0].lb == 273.15 + + assert isinstance(model.fs.unit.mole_frac_phase_comp, Var) + for j, p, i in model.fs.unit.mole_frac_phase_comp: + assert value(model.fs.unit.mole_frac_phase_comp[j, p, i]) == 0.5 + + assert_units_consistent(model) + + @pytest.mark.component + def test_structural_issues(self, model): + dt = DiagnosticsToolbox(model) + dt.assert_no_structural_warnings() + + @pytest.mark.ui + @pytest.mark.unit + def test_get_performance_contents(self, model): + perf_dict = model.fs.unit._get_performance_contents() + + assert perf_dict is None + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_initialize(self, model): + initialization_tester(model) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solve(self, model): + results = solver.solve(model) + + # Check for optimal solution + assert check_optimal_termination(results) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solution(self, model): + assert pytest.approx(101325.0, abs=1e-3) == value( + model.fs.unit.outlet.pressure[0] + ) + assert pytest.approx(368.0, abs=1e-3) == value( + model.fs.unit.outlet.temperature[0] + ) + + assert pytest.approx(1.0, abs=1e-3) == value( + model.fs.unit.control_volume.properties_out[0].flow_mol + ) + + assert pytest.approx(0.603, abs=1e-3) == value( + model.fs.unit.outlet.flow_mol_phase[0, "Liq"] + ) + assert pytest.approx(0.396, abs=1e-3) == value( + model.fs.unit.outlet.flow_mol_phase[0, "Vap"] + ) + + assert pytest.approx(0.412, abs=1e-3) == value( + model.fs.unit.outlet.mole_frac_phase_comp[0, "Liq", "benzene"] + ) + assert pytest.approx(0.588, abs=1e-3) == value( + model.fs.unit.outlet.mole_frac_phase_comp[0, "Liq", "toluene"] + ) + assert pytest.approx(0.634, abs=1e-3) == value( + model.fs.unit.outlet.mole_frac_phase_comp[0, "Vap", "benzene"] + ) + assert pytest.approx(0.366, abs=1e-3) == value( + model.fs.unit.outlet.mole_frac_phase_comp[0, "Vap", "toluene"] + ) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_numerical_issues(self, model): + dt = DiagnosticsToolbox(model) + dt.assert_no_numerical_warnings() diff --git a/idaes/models/properties/modular_properties/phase_equil/smooth_VLE.py b/idaes/models/properties/modular_properties/phase_equil/smooth_VLE.py index 463edd67e9..72c7d712cd 100644 --- a/idaes/models/properties/modular_properties/phase_equil/smooth_VLE.py +++ b/idaes/models/properties/modular_properties/phase_equil/smooth_VLE.py @@ -92,6 +92,7 @@ def phase_equil(b, phase_pair): initialize=b.temperature.value, doc="Intermediate temperature for calculating Teq", units=t_units, + bounds=b.temperature.bounds, ), ) _t1 = getattr(b, "_t1" + suffix) diff --git a/idaes/models/properties/modular_properties/state_definitions/FpTPxpc.py b/idaes/models/properties/modular_properties/state_definitions/FpTPxpc.py new file mode 100644 index 0000000000..3ba239d346 --- /dev/null +++ b/idaes/models/properties/modular_properties/state_definitions/FpTPxpc.py @@ -0,0 +1,399 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2024 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Methods for setting up FpTPxpc as the state variables in a generic property +package +""" +# TODO: Missing docstrings +# pylint: disable=missing-function-docstring + +# TODO: Look into protected access issues +# pylint: disable=protected-access + +from types import MethodType +from pyomo.environ import ( + Constraint, + Expression, + NonNegativeReals, + Var, + units as pyunits, +) +from pyomo.util.calc_var_value import calculate_variable_from_constraint + +from idaes.core import MaterialFlowBasis, MaterialBalanceType, EnergyBalanceType +from idaes.models.properties.modular_properties.base.utility import ( + get_bounds_from_config, +) + +from idaes.core.util.exceptions import ConfigurationError +import idaes.logger as idaeslog +import idaes.core.util.scaling as iscale +from .electrolyte_states import define_electrolyte_state, calculate_electrolyte_scaling + + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +def set_metadata(blk): + # This is the default assumption for state vars, so we don't need to change + # the metadata dict + pass + + +def define_state(blk): + blk.always_flash = False + + units = blk.params.get_metadata().derived_units + + # Check that only necessary state_bounds are defined + expected_keys = ["flow_mol_phase", "temperature", "pressure"] + if ( + blk.params.config.state_bounds is not None + and any(blk.params.config.state_bounds.keys()) not in expected_keys + ): + for k in blk.params.config.state_bounds.keys(): + if "mole_frac" in k: + _log.warning( + "{} - found state_bounds argument for {}." + " Mole fraction bounds are set automatically and " + "this argument will be ignored.".format(blk.name, k) + ) + elif k not in expected_keys: + raise ConfigurationError( + "{} - found unexpected state_bounds key {}. Please ensure " + "bounds are provided only for expected state variables " + "and that you have typed the variable names correctly.".format( + blk.name, k + ) + ) + + # Get bounds and initial values from config args + + f_bounds, f_init = get_bounds_from_config(blk, "flow_mol_phase", units.FLOW_MOLE) + t_bounds, t_init = get_bounds_from_config(blk, "temperature", units.TEMPERATURE) + p_bounds, p_init = get_bounds_from_config(blk, "pressure", units.PRESSURE) + + # Add state variables + blk.flow_mol_phase = Var( + blk.phase_list, + initialize=f_init, + domain=NonNegativeReals, + bounds=f_bounds, + doc="Total molar flowrate", + units=units.FLOW_MOLE, + ) + blk.mole_frac_phase_comp = Var( + blk.phase_component_set, + bounds=(1e-20, 1.001), + initialize=1 / len(blk.component_list), + doc="Mixture mole fractions", + units=pyunits.dimensionless, + ) + blk.pressure = Var( + initialize=p_init, + domain=NonNegativeReals, + bounds=p_bounds, + doc="State pressure", + units=units.PRESSURE, + ) + blk.temperature = Var( + initialize=t_init, + domain=NonNegativeReals, + bounds=t_bounds, + doc="State temperature", + units=units.TEMPERATURE, + ) + + # Add supporting variables and expressions + + @blk.Expression(doc="Total molar flow") + def flow_mol(b): + return sum(b.flow_mol_phase[p] for p in b.phase_list) + + blk.mole_frac_comp = Var( + blk.component_list, + initialize=1 / len(blk.component_list), + bounds=(1e-20, 1.001), + doc="Mole fractions", + units=pyunits.dimensionless, + ) + + @blk.Expression(blk.phase_component_set, doc="Phase-component molar flowrates") + def flow_mol_phase_comp(b, p, j): + return b.flow_mol_phase[p] * b.mole_frac_phase_comp[p, j] + + blk.phase_frac = Var( + blk.phase_list, + initialize=1.0 / len(blk.phase_list), + bounds=(0, None), + doc="Phase fractions", + units=pyunits.dimensionless, + ) + + # Add electrolyte state vars if required + if blk.params._electrolyte: + # TODO do we need special handling here? + # raise NotImplementedError() + define_electrolyte_state(blk) + + # Add supporting constraints + if blk.config.defined_state is False: + # applied at outlet only + @blk.Constraint(blk.phase_list) + def sum_mole_frac_out(b, p): + return 1.0 == sum( + blk.mole_frac_phase_comp[p, j] + for j in b.component_list + if (p, j) in b.phase_component_set + ) + + @blk.Constraint(blk.component_list, doc="Defines mole_frac_comp") + def mole_frac_comp_eq(b, j): + return b.mole_frac_comp[j] * b.flow_mol == sum( + b.mole_frac_phase_comp[p, j] * b.flow_mol_phase[p] + for p in b.phase_list + if (p, j) in b.phase_component_set + ) + + if len(blk.phase_list) == 1: + + @blk.Constraint(blk.phase_list, doc="Defines phase_frac") + def phase_fraction_constraint(b, p): + return b.phase_frac[p] == 1.0 + + else: + + def rule_phase_frac(b, p): + return b.phase_frac[p] * b.flow_mol == b.flow_mol_phase[p] + + blk.phase_fraction_constraint = Constraint(blk.phase_list, rule=rule_phase_frac) + + # ------------------------------------------------------------------------- + # General Methods + def get_material_flow_terms_FpTPxpc(b, p, j): + """Create material flow terms for control volume.""" + return b.flow_mol_phase_comp[p, j] + + blk.get_material_flow_terms = MethodType(get_material_flow_terms_FpTPxpc, blk) + + def get_enthalpy_flow_terms_FpTPxpc(b, p): + """Create enthalpy flow terms.""" + # enth_mol_phase probably does not exist when this is created + # Use try/except to build flow term if not present + try: + eflow = b._enthalpy_flow_term + except AttributeError: + + def rule_eflow(b, p): + return b.flow_mol_phase[p] * b.enth_mol_phase[p] + + eflow = b._enthalpy_flow_term = Expression(b.phase_list, rule=rule_eflow) + return eflow[p] + + blk.get_enthalpy_flow_terms = MethodType(get_enthalpy_flow_terms_FpTPxpc, blk) + + def get_material_density_terms_FpTPxpc(b, p, j): + """Create material density terms.""" + # dens_mol_phase probably does not exist when this is created + # Use try/except to build term if not present + try: + mdens = b._material_density_term + except AttributeError: + + def rule_mdens(b, p, j): + return b.dens_mol_phase[p] * b.mole_frac_phase_comp[p, j] + + mdens = b._material_density_term = Expression( + b.phase_component_set, rule=rule_mdens + ) + return mdens[p, j] + + blk.get_material_density_terms = MethodType(get_material_density_terms_FpTPxpc, blk) + + def get_energy_density_terms_FpTPxpc(b, p): + """Create energy density terms.""" + # Density and energy terms probably do not exist when this is created + # Use try/except to build term if not present + try: + edens = b._energy_density_term + except AttributeError: + + def rule_edens(b, p): + return b.dens_mol_phase[p] * b.energy_internal_mol_phase[p] + + edens = b._energy_density_term = Expression(b.phase_list, rule=rule_edens) + return edens[p] + + blk.get_energy_density_terms = MethodType(get_energy_density_terms_FpTPxpc, blk) + + def default_material_balance_type_FpTPxpc(): + return MaterialBalanceType.componentTotal + + blk.default_material_balance_type = default_material_balance_type_FpTPxpc + + def default_energy_balance_type_FpTPxpc(): + return EnergyBalanceType.enthalpyTotal + + blk.default_energy_balance_type = default_energy_balance_type_FpTPxpc + + def get_material_flow_basis_FpTPxpc(): + return MaterialFlowBasis.molar + + blk.get_material_flow_basis = get_material_flow_basis_FpTPxpc + + def define_state_vars_FpTPxpc(b): + """Define state vars.""" + return { + "flow_mol_phase": b.flow_mol_phase, + "mole_frac_phase_comp": b.mole_frac_phase_comp, + "temperature": b.temperature, + "pressure": b.pressure, + } + + blk.define_state_vars = MethodType(define_state_vars_FpTPxpc, blk) + + def define_display_vars_FpTPxpc(b): + """Define display vars.""" + return { + "Total Molar Flowrate": b.flow_mol_phase, + "Total Mole Fraction": b.mole_frac_phase_comp, + "Temperature": b.temperature, + "Pressure": b.pressure, + } + + blk.define_display_vars = MethodType(define_display_vars_FpTPxpc, blk) + + +# TODO flash initialization---needs to be a separate method because +# we need to build a ControlVolume0D to enforce material balances +def state_initialization(b): + """ + Initialize state variables that need to be initialize irrespective + of whether or not there's a flash. + """ + for j in b.component_list: + # Linear in mole_frac_comp when all other variables + # are fixed---this should converge in one iteration + calculate_variable_from_constraint(b.mole_frac_comp[j], b.mole_frac_comp_eq[j]) + + +def define_default_scaling_factors(b): + """ + Method to set default scaling factors for the property package. Scaling + factors are based on the default initial value for each variable provided + in the state_bounds config argument. + """ + # Get bounds and initial values from config args + units = b.get_metadata().derived_units + state_bounds = b.config.state_bounds + + if state_bounds is None: + return + + try: + f_bounds = state_bounds["flow_mol_phase"] + if len(f_bounds) == 4: + f_init = pyunits.convert_value( + f_bounds[1], from_units=f_bounds[3], to_units=units.FLOW_MOLE + ) + else: + f_init = f_bounds[1] + except KeyError: + f_init = 1 + + try: + p_bounds = state_bounds["pressure"] + if len(p_bounds) == 4: + p_init = pyunits.convert_value( + p_bounds[1], from_units=p_bounds[3], to_units=units.PRESSURE + ) + else: + p_init = p_bounds[1] + except KeyError: + p_init = 1 + + try: + t_bounds = state_bounds["temperature"] + if len(t_bounds) == 4: + t_init = pyunits.convert_value( + t_bounds[1], from_units=t_bounds[3], to_units=units.TEMPERATURE + ) + else: + t_init = t_bounds[1] + except KeyError: + t_init = 1 + + # Set default scaling factors + b.set_default_scaling("flow_mol", 1 / f_init) + b.set_default_scaling("flow_mol_phase", 1 / f_init) + b.set_default_scaling("flow_mol_comp", 1 / f_init) + b.set_default_scaling("flow_mol_phase_comp", 1 / f_init) + b.set_default_scaling("pressure", 1 / p_init) + b.set_default_scaling("temperature", 1 / t_init) + + +def calculate_scaling_factors(blk): + sf_flow_phase = {} + for p, v in blk.flow_mol_phase.items(): + sf_flow_phase[p] = iscale.get_scaling_factor( + blk.flow_mol_phase[p], default=1, warning=True + ) + sf_mf = {} + for i, v in blk.mole_frac_phase_comp.items(): + sf_mf[i] = iscale.get_scaling_factor(v, default=1e3, warning=True) + + for i in blk.mole_frac_comp: + sf = iscale.get_scaling_factor(blk.mole_frac_comp[i]) + if sf is None: + sf = min(sf_mf[p, j] for j, p in blk.phase_component_set if i == j) + iscale.set_scaling_factor(sf, blk.mole_frac_comp[i]) + iscale.constraint_scaling_transform( + blk.mole_frac_comp_eq[i], sf, overwrite=False + ) + + if blk.config.defined_state is False: + for p in blk.phase_list: + sf_eqn = min(sf_mf[p, j] for j in blk.components_in_phase(p)) + iscale.constraint_scaling_transform( + blk.sum_mole_frac_out[p], sf_eqn, overwrite=False + ) + + if len(blk.phase_list) == 1: + for p in blk.phase_list: + iscale.constraint_scaling_transform( + blk.phase_fraction_constraint[p], 1, overwrite=False + ) + else: + for p in blk.phase_list: + iscale.constraint_scaling_transform( + blk.phase_fraction_constraint[p], sf_flow_phase[p], overwrite=False + ) + + if blk.params._electrolyte: + # raise NotImplementedError() + calculate_electrolyte_scaling(blk) + + +do_not_initialize = ["sum_mole_frac_out"] + + +class FpTPxpc(object): + """Total flow, temperature, pressure, mole fraction state.""" + + set_metadata = set_metadata + define_state = define_state + state_initialization = state_initialization + do_not_initialize = do_not_initialize + define_default_scaling_factors = define_default_scaling_factors + calculate_scaling_factors = calculate_scaling_factors diff --git a/idaes/models/properties/modular_properties/state_definitions/tests/test_FpTPxpc.py b/idaes/models/properties/modular_properties/state_definitions/tests/test_FpTPxpc.py new file mode 100644 index 0000000000..43691c8ac7 --- /dev/null +++ b/idaes/models/properties/modular_properties/state_definitions/tests/test_FpTPxpc.py @@ -0,0 +1,1556 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2024 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Tests for FpTPxpc state formulation + +Authors: Tanner Polley +""" + +import pytest +import re +from sys import modules + +from pyomo.environ import ConcreteModel, Constraint, Expression, Var, units as pyunits +from pyomo.util.check_units import check_units_equivalent, assert_units_consistent + +# Need define_default_scaling_factors, even though it is not used directly +from idaes.models.properties.modular_properties.state_definitions.FpTPxpc import ( + FpTPxpc, + define_state, + set_metadata, +) +from idaes.core import ( + FlowsheetBlock, + MaterialFlowBasis, + MaterialBalanceType, + EnergyBalanceType, + declare_process_block_class, + Solvent, + Solute, + AqueousPhase, + VaporPhase, +) +from idaes.models.properties.modular_properties.base.generic_property import ( + GenericParameterData, + GenericParameterBlock, +) +from idaes.models.properties.modular_properties.base.tests.dummy_eos import DummyEoS +from idaes.models.properties.modular_properties.pure.Perrys import Perrys +from idaes.models.properties.modular_properties.pure.NIST import NIST +from idaes.models.properties.modular_properties.phase_equil.forms import fugacity +from idaes.models.properties.modular_properties.eos.ideal import Ideal +from idaes.models.properties.modular_properties.phase_equil import SmoothVLE +from idaes.models.properties.modular_properties.phase_equil.bubble_dew import ( + IdealBubbleDew, +) +from idaes.core.util.exceptions import ConfigurationError +import idaes.logger as idaeslog +from idaes.core.util.model_statistics import degrees_of_freedom, large_residuals_set + + +# Note: The state definition is set by importing functions for the relevant module above +@declare_process_block_class("DummyParameterBlock") +class DummyParameterData(GenericParameterData): + pass + + +@pytest.mark.unit +def test_set_metadata(): + assert set_metadata(None) is None + + +class TestInvalidBounds(object): + @pytest.mark.unit + def test_bad_name(self): + m = ConcreteModel() + + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={"p1": {"equation_of_state": DummyEoS}}, + state_definition=modules[__name__], + pressure_ref=100000.0, + temperature_ref=300, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + state_bounds={"foo": (None, None, None)}, + ) + + with pytest.raises( + ConfigurationError, + match=re.escape( + "props[1] - found unexpected state_bounds key foo. " + "Please ensure bounds are provided only for expected state " + "variables and that you have typed the variable names " + "correctly." + ), + ): + # Build state block + m.props = m.params.build_state_block([1], defined_state=True) + + @pytest.mark.unit + def test_mole_frac(self, caplog): + m = ConcreteModel() + + caplog.set_level( + idaeslog.WARNING, logger=("idaes.models.properties.modular_properties.") + ) + + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={"p1": {"equation_of_state": DummyEoS}}, + state_definition=modules[__name__], + pressure_ref=100000.0, + temperature_ref=300, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + state_bounds={"mole_frac_phase_comp": (None, None, None)}, + ) + + # Create a dummy state block + m.props = m.params.build_state_block([1], defined_state=True) + + assert ( + "props[1] - found state_bounds argument for mole_frac_phase_comp." + " Mole fraction bounds are set automatically and " + "this argument will be ignored." in caplog.text + ) + + +class Test1PhaseDefinedStateFalseNoBounds(object): + # Test define_state method with no bounds and defined_State = False + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={"p1": {"equation_of_state": DummyEoS}}, + state_definition=modules[__name__], + pressure_ref=100000.0, + temperature_ref=300, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + ) + + # Build state block + m.props = m.params.build_state_block([1], defined_state=False) + + # Add necessary variables that would be built by other methods + m.props[1].dens_mol_phase = Var(m.params.phase_list, initialize=1) + m.props[1].enth_mol_phase = Var(m.params.phase_list, initialize=1) + + return m + + @pytest.mark.unit + def test_always_flash(self, frame): + define_state(frame.props[1]) + + assert not frame.props[1].always_flash + + @pytest.mark.unit + def test_vars(self, frame): + # Check that all necessary variables have been constructed and have + # the correct values + + # flow_mol_phase Var line 97 + assert isinstance(frame.props[1].flow_mol_phase, Var) + assert len(frame.props[1].flow_mol_phase) == 1 + for p in frame.props[1].flow_mol_phase: + assert p in frame.params.phase_list + assert frame.props[1].flow_mol_phase[p].value is None + assert check_units_equivalent( + frame.props[1].flow_mol_phase, pyunits.mol / pyunits.s + ) + + # mole_frac_phase_comp Var line 105 + assert isinstance(frame.props[1].mole_frac_phase_comp, Var) + assert len(frame.props[1].mole_frac_phase_comp) == 3 + for p, i in frame.props[1].mole_frac_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert frame.props[1].mole_frac_phase_comp[p, i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_phase_comp, None) + + # pressure Var line 112 + assert isinstance(frame.props[1].pressure, Var) + assert frame.props[1].pressure.value is None + assert check_units_equivalent(frame.props[1].pressure, pyunits.Pa) + + # temperature Var line 119 + assert isinstance(frame.props[1].temperature, Var) + assert frame.props[1].temperature.value is None + assert check_units_equivalent(frame.props[1].temperature, pyunits.K) + + # mole_frac_comp Var line 158 + assert isinstance(frame.props[1].mole_frac_comp, Var) + assert len(frame.props[1].mole_frac_comp) == 3 + for i in frame.props[1].mole_frac_comp: + assert i in frame.props[1].params.component_list + assert frame.props[1].mole_frac_comp[i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_comp, None) + + # phase_frac Expression line 168 + assert isinstance(frame.props[1].phase_frac, Var) + assert len(frame.props[1].phase_frac) == 1 + for i in frame.props[1].phase_frac: + assert i in frame.props[1].params.phase_list + assert frame.props[1].phase_frac[i].value == 1 + assert check_units_equivalent(frame.props[1].phase_frac, None) + + @pytest.mark.unit + def test_expressions(self, frame): + + # flow_mol Expression line 129 + assert isinstance(frame.props[1].flow_mol, Expression) + assert len(frame.props[1].flow_mol) == 1 + assert str(frame.props[1].flow_mol.expr) == str( + sum( + frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + ) + ) + + # flow_mol_phase_comp Expression line 133 + assert isinstance(frame.props[1].flow_mol_phase_comp, Expression) + assert len(frame.props[1].flow_mol_phase_comp) == 3 + for p, i in frame.props[1].flow_mol_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert str(frame.props[1].flow_mol_phase_comp[p, i].expr) == str( + frame.props[1].flow_mol_phase[p] + * frame.props[1].mole_frac_phase_comp[p, i] + ) + + @pytest.mark.unit + def test_constraints(self, frame): + # Check that the correct constraints are present + + # sum_mole_frac_out Constraint line 179 + assert isinstance(frame.props[1].sum_mole_frac_out, Constraint) + assert len(frame.props[1].sum_mole_frac_out) == 1 + for p in frame.props[1].phase_frac: + assert p in frame.params.phase_list + assert str(frame.props[1].sum_mole_frac_out[p].body) == str( + sum( + frame.props[1].mole_frac_phase_comp[p, i] + for i in frame.props[1].params.component_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # mole_frac_comp_eq Constraint line 183 + assert isinstance(frame.props[1].mole_frac_comp_eq, Constraint) + assert len(frame.props[1].mole_frac_comp_eq) == 3 + for i in frame.props[1].mole_frac_comp_eq: + assert str(frame.props[1].mole_frac_comp_eq[i].body) == str( + frame.props[1].mole_frac_comp[i] * frame.props[1].flow_mol + - sum( + frame.props[1].mole_frac_phase_comp[p, i] + * frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # phase_fraction_constraint Constraint line 200 + assert isinstance(frame.props[1].phase_fraction_constraint, Constraint) + assert len(frame.props[1].phase_fraction_constraint) == 1 + for p in frame.props[1].phase_fraction_constraint: + assert p in frame.params.phase_list + assert str(frame.props[1].phase_fraction_constraint[p].body) == str( + frame.props[1].phase_frac[p] + ) + + assert_units_consistent(frame.props[1]) + + +class Test1PhaseDefinedStateTrueWithBounds(object): + # Test define_state method with bounds and defined_State = True + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + # Create a dummy parameter block + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={"p1": {"equation_of_state": DummyEoS}}, + state_definition=modules[__name__], + pressure_ref=100000.0, + temperature_ref=300, + state_bounds={ + "flow_mol_phase": (0, 100, 200), + "temperature": (290, 345, 400), + "pressure": (100000.0, 300000.0, 500000.0), + }, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + ) + + # Build state block + m.props = m.params.build_state_block([1], defined_state=True, has_phase_equilibrium=False) + + # Add necessary variables that would be built by other methods + m.props[1].dens_mol_phase = Var(m.params.phase_list, initialize=1) + m.props[1].enth_mol_phase = Var(m.params.phase_list, initialize=1) + + return m + + @pytest.mark.unit + def test_always_flash(self, frame): + define_state(frame.props[1]) + + assert not frame.props[1].always_flash + + @pytest.mark.unit + def test_vars(self, frame): + # Check that all necessary variables have been constructed and have + # the correct values + + # flow_mol_phase Var line 97 + assert isinstance(frame.props[1].flow_mol_phase, Var) + assert len(frame.props[1].flow_mol_phase) == 1 + for p in frame.props[1].flow_mol_phase: + assert p in frame.params.phase_list + assert frame.props[1].flow_mol_phase[p].value == 100 + assert frame.props[1].flow_mol_phase[p].lb == 0 + assert frame.props[1].flow_mol_phase[p].ub == 200 + assert check_units_equivalent( + frame.props[1].flow_mol_phase, pyunits.mol / pyunits.s + ) + + # mole_frac_phase_comp Var line 105 + assert isinstance(frame.props[1].mole_frac_phase_comp, Var) + assert len(frame.props[1].mole_frac_phase_comp) == 3 + for p, i in frame.props[1].mole_frac_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert frame.props[1].mole_frac_phase_comp[p, i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_phase_comp, None) + + # pressure Var line 112 + assert isinstance(frame.props[1].pressure, Var) + assert frame.props[1].pressure.value == 3e5 + assert frame.props[1].pressure.lb == 1e5 + assert frame.props[1].pressure.ub == 5e5 + assert check_units_equivalent(frame.props[1].pressure, pyunits.Pa) + + # temperature Var line 119 + assert isinstance(frame.props[1].temperature, Var) + assert frame.props[1].temperature.value == 345 + assert frame.props[1].temperature.lb == 290 + assert frame.props[1].temperature.ub == 400 + assert check_units_equivalent(frame.props[1].temperature, pyunits.K) + + # mole_frac_comp Var line 158 + assert isinstance(frame.props[1].mole_frac_comp, Var) + assert len(frame.props[1].mole_frac_comp) == 3 + for i in frame.props[1].mole_frac_comp: + assert i in frame.props[1].params.component_list + assert frame.props[1].mole_frac_comp[i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_comp, None) + + # phase_frac Expression line 168 + assert isinstance(frame.props[1].phase_frac, Var) + assert len(frame.props[1].phase_frac) == 1 + for i in frame.props[1].phase_frac: + assert i in frame.props[1].params.phase_list + assert frame.props[1].phase_frac[i].value == 1 + assert check_units_equivalent(frame.props[1].phase_frac, None) + + @pytest.mark.unit + def test_expressions(self, frame): + + # flow_mol Expression line 129 + assert isinstance(frame.props[1].flow_mol, Expression) + assert len(frame.props[1].flow_mol) == 1 + assert str(frame.props[1].flow_mol.expr) == str( + sum( + frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + ) + ) + + # flow_mol_phase_comp Expression line 133 + assert isinstance(frame.props[1].flow_mol_phase_comp, Expression) + assert len(frame.props[1].flow_mol_phase_comp) == 3 + for p, i in frame.props[1].flow_mol_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert str(frame.props[1].flow_mol_phase_comp[p, i].expr) == str( + frame.props[1].flow_mol_phase[p] + * frame.props[1].mole_frac_phase_comp[p, i] + ) + + @pytest.mark.unit + def test_constraints(self, frame): + # Check that the correct constraints are present + + # mole_frac_comp_eq Constraint line 183 + assert isinstance(frame.props[1].mole_frac_comp_eq, Constraint) + assert len(frame.props[1].mole_frac_comp_eq) == 3 + for i in frame.props[1].mole_frac_comp_eq: + assert str(frame.props[1].mole_frac_comp_eq[i].body) == str( + frame.props[1].mole_frac_comp[i] * frame.props[1].flow_mol + - sum( + frame.props[1].mole_frac_phase_comp[p, i] + * frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # phase_fraction_constraint Constraint line 200 + assert isinstance(frame.props[1].phase_fraction_constraint, Constraint) + assert len(frame.props[1].phase_fraction_constraint) == 1 + for p in frame.props[1].phase_fraction_constraint: + assert p in frame.params.phase_list + assert str(frame.props[1].phase_fraction_constraint[p].body) == str( + frame.props[1].phase_frac[p] + ) + + assert_units_consistent(frame.props[1]) + + +class Test2PhaseDefinedStateFalseNoBounds(object): + # Test define_state method with no bounds and defined_State = False + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + # Create a dummy parameter block + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={ + "p1": {"equation_of_state": DummyEoS}, + "p2": {"equation_of_state": DummyEoS}, + }, + state_definition=modules[__name__], + pressure_ref=100000.0, + temperature_ref=300, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + ) + + # Build state block + m.props = m.params.build_state_block([1], defined_state=False) + + # Add necessary variables that would be built by other methods + m.props[1].dens_mol_phase = Var(m.params.phase_list, initialize=1) + m.props[1].enth_mol_phase = Var(m.params.phase_list, initialize=1) + + return m + + @pytest.mark.unit + def test_always_flash(self, frame): + define_state(frame.props[1]) + + assert not frame.props[1].always_flash + + @pytest.mark.unit + def test_vars(self, frame): + # Check that all necessary variables have been constructed and have + # the correct values + + # flow_mol_phase Var line 97 + assert isinstance(frame.props[1].flow_mol_phase, Var) + assert len(frame.props[1].flow_mol_phase) == 2 + for p in frame.props[1].flow_mol_phase: + assert p in frame.params.phase_list + assert frame.props[1].flow_mol_phase[p].value is None + assert check_units_equivalent( + frame.props[1].flow_mol_phase, pyunits.mol / pyunits.s + ) + + # mole_frac_phase_comp Var line 105 + assert isinstance(frame.props[1].mole_frac_phase_comp, Var) + assert len(frame.props[1].mole_frac_phase_comp) == 6 + for p, i in frame.props[1].mole_frac_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert frame.props[1].mole_frac_phase_comp[p, i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_phase_comp, None) + + # pressure Var line 112 + assert isinstance(frame.props[1].pressure, Var) + assert frame.props[1].pressure.value is None + assert check_units_equivalent(frame.props[1].pressure, pyunits.Pa) + + # temperature Var line 119 + assert isinstance(frame.props[1].temperature, Var) + assert frame.props[1].temperature.value is None + assert check_units_equivalent(frame.props[1].temperature, pyunits.K) + + # mole_frac_comp Var line 158 + assert isinstance(frame.props[1].mole_frac_comp, Var) + assert len(frame.props[1].mole_frac_comp) == 3 + for i in frame.props[1].mole_frac_comp: + assert i in frame.props[1].params.component_list + assert frame.props[1].mole_frac_comp[i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_comp, None) + + # phase_frac Expression line 168 + assert isinstance(frame.props[1].phase_frac, Var) + assert len(frame.props[1].phase_frac) == 2 + for i in frame.props[1].phase_frac: + assert i in frame.props[1].params.phase_list + assert frame.props[1].phase_frac[i].value == 1 / 2 + assert check_units_equivalent(frame.props[1].phase_frac, None) + + @pytest.mark.unit + def test_expressions(self, frame): + + # flow_mol Expression line 129 + assert isinstance(frame.props[1].flow_mol, Expression) + assert len(frame.props[1].flow_mol) == 1 + assert str(frame.props[1].flow_mol.expr) == str( + sum( + frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + ) + ) + + # flow_mol_phase_comp Expression line 133 + assert isinstance(frame.props[1].flow_mol_phase_comp, Expression) + assert len(frame.props[1].flow_mol_phase_comp) == 6 + for p, i in frame.props[1].flow_mol_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert str(frame.props[1].flow_mol_phase_comp[p, i].expr) == str( + frame.props[1].flow_mol_phase[p] + * frame.props[1].mole_frac_phase_comp[p, i] + ) + + @pytest.mark.unit + def test_constraints(self, frame): + # Check that the correct constraints are present + + # sum_mole_frac_out Constraint line 179 + assert isinstance(frame.props[1].sum_mole_frac_out, Constraint) + assert len(frame.props[1].sum_mole_frac_out) == 2 + for p in frame.props[1].phase_frac: + assert p in frame.params.phase_list + assert str(frame.props[1].sum_mole_frac_out[p].body) == str( + sum( + frame.props[1].mole_frac_phase_comp[p, i] + for i in frame.props[1].params.component_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # mole_frac_comp_eq Constraint line 183 + assert isinstance(frame.props[1].mole_frac_comp_eq, Constraint) + assert len(frame.props[1].mole_frac_comp_eq) == 3 + for i in frame.props[1].mole_frac_comp_eq: + assert str(frame.props[1].mole_frac_comp_eq[i].body) == str( + frame.props[1].mole_frac_comp[i] * frame.props[1].flow_mol + - sum( + frame.props[1].mole_frac_phase_comp[p, i] + * frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # phase_fraction_constraint Constraint line 200 + assert isinstance(frame.props[1].phase_fraction_constraint, Constraint) + assert len(frame.props[1].phase_fraction_constraint) == 2 + for p in frame.props[1].phase_fraction_constraint: + assert p in frame.params.phase_list + assert str(frame.props[1].phase_fraction_constraint[p].body) == str( + frame.props[1].phase_frac[p] * frame.props[1].flow_mol + - frame.props[1].flow_mol_phase[p] + ) + + assert_units_consistent(frame.props[1]) + + +class Test2PhaseDefinedStateTrueWithBounds(object): + # Test define_state method with no bounds and defined_State = False + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + # Create a dummy parameter block + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={ + "p1": {"equation_of_state": DummyEoS}, + "p2": {"equation_of_state": DummyEoS}, + }, + state_definition=modules[__name__], + pressure_ref=100000.0, + temperature_ref=300, + state_bounds={ + "flow_mol_phase": (0, 100, 200), + "temperature": (290, 345, 400), + "pressure": (100000.0, 300000.0, 500000.0), + }, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + ) + + # Build state block + m.props = m.params.build_state_block([1], defined_state=True) + + # Add necessary variables that would be built by other methods + m.props[1].dens_mol_phase = Var(m.params.phase_list, initialize=1) + m.props[1].enth_mol_phase = Var(m.params.phase_list, initialize=1) + + return m + + @pytest.mark.unit + def test_always_flash(self, frame): + define_state(frame.props[1]) + + assert not frame.props[1].always_flash + + @pytest.mark.unit + def test_vars(self, frame): + # Check that all necessary variables have been constructed and have + # the correct values + + # flow_mol_phase Var line 97 + assert isinstance(frame.props[1].flow_mol_phase, Var) + assert len(frame.props[1].flow_mol_phase) == 2 + for p in frame.props[1].flow_mol_phase: + assert p in frame.params.phase_list + assert frame.props[1].flow_mol_phase[p].value == 100 + assert frame.props[1].flow_mol_phase[p].lb == 0 + assert frame.props[1].flow_mol_phase[p].ub == 200 + assert check_units_equivalent( + frame.props[1].flow_mol_phase, pyunits.mol / pyunits.s + ) + + # mole_frac_phase_comp Var line 105 + assert isinstance(frame.props[1].mole_frac_phase_comp, Var) + assert len(frame.props[1].mole_frac_phase_comp) == 6 + for p, i in frame.props[1].mole_frac_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert frame.props[1].mole_frac_phase_comp[p, i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_phase_comp, None) + + # pressure Var line 112 + assert isinstance(frame.props[1].pressure, Var) + assert frame.props[1].pressure.value == 3e5 + assert frame.props[1].pressure.lb == 1e5 + assert frame.props[1].pressure.ub == 5e5 + assert check_units_equivalent(frame.props[1].pressure, pyunits.Pa) + + # temperature Var line 119 + assert isinstance(frame.props[1].temperature, Var) + assert frame.props[1].temperature.value == 345 + assert frame.props[1].temperature.lb == 290 + assert frame.props[1].temperature.ub == 400 + assert check_units_equivalent(frame.props[1].temperature, pyunits.K) + + # mole_frac_comp Var line 158 + assert isinstance(frame.props[1].mole_frac_comp, Var) + assert len(frame.props[1].mole_frac_comp) == 3 + for i in frame.props[1].mole_frac_comp: + assert i in frame.props[1].params.component_list + assert frame.props[1].mole_frac_comp[i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_comp, None) + + # phase_frac Expression line 168 + assert isinstance(frame.props[1].phase_frac, Var) + assert len(frame.props[1].phase_frac) == 2 + for i in frame.props[1].phase_frac: + assert i in frame.props[1].params.phase_list + assert frame.props[1].phase_frac[i].value == 1 / 2 + assert check_units_equivalent(frame.props[1].phase_frac, None) + + @pytest.mark.unit + def test_expressions(self, frame): + + # flow_mol Expression line 129 + assert isinstance(frame.props[1].flow_mol, Expression) + assert len(frame.props[1].flow_mol) == 1 + assert str(frame.props[1].flow_mol.expr) == str( + sum( + frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + ) + ) + + # flow_mol_phase_comp Expression line 133 + assert isinstance(frame.props[1].flow_mol_phase_comp, Expression) + assert len(frame.props[1].flow_mol_phase_comp) == 6 + for p, i in frame.props[1].flow_mol_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert str(frame.props[1].flow_mol_phase_comp[p, i].expr) == str( + frame.props[1].flow_mol_phase[p] + * frame.props[1].mole_frac_phase_comp[p, i] + ) + + @pytest.mark.unit + def test_constraints(self, frame): + # Check that the correct constraints are present + + # mole_frac_comp_eq Constraint line 183 + assert isinstance(frame.props[1].mole_frac_comp_eq, Constraint) + assert len(frame.props[1].mole_frac_comp_eq) == 3 + for i in frame.props[1].mole_frac_comp_eq: + assert str(frame.props[1].mole_frac_comp_eq[i].body) == str( + frame.props[1].mole_frac_comp[i] * frame.props[1].flow_mol + - sum( + frame.props[1].mole_frac_phase_comp[p, i] + * frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # phase_fraction_constraint Constraint line 200 + assert isinstance(frame.props[1].phase_fraction_constraint, Constraint) + assert len(frame.props[1].phase_fraction_constraint) == 2 + for p in frame.props[1].phase_fraction_constraint: + assert p in frame.params.phase_list + assert str(frame.props[1].phase_fraction_constraint[p].body) == str( + frame.props[1].phase_frac[p] * frame.props[1].flow_mol + - frame.props[1].flow_mol_phase[p] + ) + + assert_units_consistent(frame.props[1]) + + +class Test3PhaseDefinedStateFalseNoBounds(object): + # Test define_state method with no bounds and defined_State = False + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + # Create a dummy parameter block + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={ + "p1": {"equation_of_state": DummyEoS}, + "p2": {"equation_of_state": DummyEoS}, + "p3": {"equation_of_state": DummyEoS}, + }, + state_definition=modules[__name__], + pressure_ref=100000.0, + temperature_ref=300, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + ) + + # Build state block + m.props = m.params.build_state_block([1], defined_state=False) + + # Add necessary variables that would be built by other methods + m.props[1].dens_mol_phase = Var(m.params.phase_list, initialize=1) + m.props[1].enth_mol_phase = Var(m.params.phase_list, initialize=1) + + return m + + @pytest.mark.unit + def test_always_flash(self, frame): + define_state(frame.props[1]) + + assert not frame.props[1].always_flash + + @pytest.mark.unit + def test_vars(self, frame): + # Check that all necessary variables have been constructed and have + # the correct values + + # flow_mol_phase Var line 97 + assert isinstance(frame.props[1].flow_mol_phase, Var) + assert len(frame.props[1].flow_mol_phase) == 3 + for p in frame.props[1].flow_mol_phase: + assert p in frame.params.phase_list + assert frame.props[1].flow_mol_phase[p].value is None + assert check_units_equivalent( + frame.props[1].flow_mol_phase, pyunits.mol / pyunits.s + ) + + # mole_frac_phase_comp Var line 105 + assert isinstance(frame.props[1].mole_frac_phase_comp, Var) + assert len(frame.props[1].mole_frac_phase_comp) == 9 + for p, i in frame.props[1].mole_frac_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert frame.props[1].mole_frac_phase_comp[p, i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_phase_comp, None) + + # pressure Var line 112 + assert isinstance(frame.props[1].pressure, Var) + assert frame.props[1].pressure.value is None + assert check_units_equivalent(frame.props[1].pressure, pyunits.Pa) + + # temperature Var line 119 + assert isinstance(frame.props[1].temperature, Var) + assert frame.props[1].temperature.value is None + assert check_units_equivalent(frame.props[1].temperature, pyunits.K) + + # mole_frac_comp Var line 158 + assert isinstance(frame.props[1].mole_frac_comp, Var) + assert len(frame.props[1].mole_frac_comp) == 3 + for i in frame.props[1].mole_frac_comp: + assert i in frame.props[1].params.component_list + assert frame.props[1].mole_frac_comp[i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_comp, None) + + # phase_frac Expression line 168 + assert isinstance(frame.props[1].phase_frac, Var) + assert len(frame.props[1].phase_frac) == 3 + for i in frame.props[1].phase_frac: + assert i in frame.props[1].params.phase_list + assert frame.props[1].phase_frac[i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].phase_frac, None) + + @pytest.mark.unit + def test_expressions(self, frame): + + # flow_mol Expression line 129 + assert isinstance(frame.props[1].flow_mol, Expression) + assert len(frame.props[1].flow_mol) == 1 + assert str(frame.props[1].flow_mol.expr) == str( + sum( + frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + ) + ) + + # flow_mol_phase_comp Expression line 133 + assert isinstance(frame.props[1].flow_mol_phase_comp, Expression) + assert len(frame.props[1].flow_mol_phase_comp) == 9 + for p, i in frame.props[1].flow_mol_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert str(frame.props[1].flow_mol_phase_comp[p, i].expr) == str( + frame.props[1].flow_mol_phase[p] + * frame.props[1].mole_frac_phase_comp[p, i] + ) + + @pytest.mark.unit + def test_constraints(self, frame): + # Check that the correct constraints are present + + # sum_mole_frac_out Constraint line 179 + assert isinstance(frame.props[1].sum_mole_frac_out, Constraint) + assert len(frame.props[1].sum_mole_frac_out) == 3 + for p in frame.props[1].phase_frac: + assert p in frame.params.phase_list + assert str(frame.props[1].sum_mole_frac_out[p].body) == str( + sum( + frame.props[1].mole_frac_phase_comp[p, i] + for i in frame.props[1].params.component_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # mole_frac_comp_eq Constraint line 183 + assert isinstance(frame.props[1].mole_frac_comp_eq, Constraint) + assert len(frame.props[1].mole_frac_comp_eq) == 3 + for i in frame.props[1].mole_frac_comp_eq: + assert str(frame.props[1].mole_frac_comp_eq[i].body) == str( + frame.props[1].mole_frac_comp[i] * frame.props[1].flow_mol + - sum( + frame.props[1].mole_frac_phase_comp[p, i] + * frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # phase_fraction_constraint Constraint line 200 + assert isinstance(frame.props[1].phase_fraction_constraint, Constraint) + assert len(frame.props[1].phase_fraction_constraint) == 3 + for p in frame.props[1].phase_fraction_constraint: + assert p in frame.params.phase_list + assert str(frame.props[1].phase_fraction_constraint[p].body) == str( + frame.props[1].phase_frac[p] * frame.props[1].flow_mol + - frame.props[1].flow_mol_phase[p] + ) + + assert_units_consistent(frame.props[1]) + + +class Test3PhaseDefinedStateTrueWithBounds(object): + # Test define_state method with no bounds and defined_State = False + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + # Create a dummy parameter block + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={ + "p1": {"equation_of_state": DummyEoS}, + "p2": {"equation_of_state": DummyEoS}, + "p3": {"equation_of_state": DummyEoS}, + }, + state_definition=modules[__name__], + pressure_ref=100000.0, + temperature_ref=300, + state_bounds={ + "flow_mol_phase": (0, 100, 200), + "temperature": (290, 345, 400), + "pressure": (100000.0, 300000.0, 500000.0), + }, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + ) + + # Build state block + m.props = m.params.build_state_block([1], defined_state=True) + + # Add necessary variables that would be built by other methods + m.props[1].dens_mol_phase = Var(m.params.phase_list, initialize=1) + m.props[1].enth_mol_phase = Var(m.params.phase_list, initialize=1) + + return m + + @pytest.mark.unit + def test_always_flash(self, frame): + define_state(frame.props[1]) + + assert not frame.props[1].always_flash + + @pytest.mark.unit + def test_vars(self, frame): + # Check that all necessary variables have been constructed and have + # the correct values + + # flow_mol_phase Var line 97 + assert isinstance(frame.props[1].flow_mol_phase, Var) + assert len(frame.props[1].flow_mol_phase) == 3 + for p in frame.props[1].flow_mol_phase: + assert p in frame.params.phase_list + assert frame.props[1].flow_mol_phase[p].value == 100 + assert frame.props[1].flow_mol_phase[p].lb == 0 + assert frame.props[1].flow_mol_phase[p].ub == 200 + assert check_units_equivalent( + frame.props[1].flow_mol_phase, pyunits.mol / pyunits.s + ) + + # mole_frac_phase_comp Var line 105 + assert isinstance(frame.props[1].mole_frac_phase_comp, Var) + assert len(frame.props[1].mole_frac_phase_comp) == 9 + for p, i in frame.props[1].mole_frac_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert frame.props[1].mole_frac_phase_comp[p, i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_phase_comp, None) + + # pressure Var line 112 + assert isinstance(frame.props[1].pressure, Var) + assert frame.props[1].pressure.value == 3e5 + assert frame.props[1].pressure.lb == 1e5 + assert frame.props[1].pressure.ub == 5e5 + assert check_units_equivalent(frame.props[1].pressure, pyunits.Pa) + + # temperature Var line 119 + assert isinstance(frame.props[1].temperature, Var) + assert frame.props[1].temperature.value == 345 + assert frame.props[1].temperature.lb == 290 + assert frame.props[1].temperature.ub == 400 + assert check_units_equivalent(frame.props[1].temperature, pyunits.K) + + # mole_frac_comp Var line 158 + assert isinstance(frame.props[1].mole_frac_comp, Var) + assert len(frame.props[1].mole_frac_comp) == 3 + for i in frame.props[1].mole_frac_comp: + assert i in frame.props[1].params.component_list + assert frame.props[1].mole_frac_comp[i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].mole_frac_comp, None) + + # phase_frac Expression line 168 + assert isinstance(frame.props[1].phase_frac, Var) + assert len(frame.props[1].phase_frac) == 3 + for i in frame.props[1].phase_frac: + assert i in frame.props[1].params.phase_list + assert frame.props[1].phase_frac[i].value == 1 / 3 + assert check_units_equivalent(frame.props[1].phase_frac, None) + + @pytest.mark.unit + def test_expressions(self, frame): + + # flow_mol Expression line 129 + assert isinstance(frame.props[1].flow_mol, Expression) + assert len(frame.props[1].flow_mol) == 1 + assert str(frame.props[1].flow_mol.expr) == str( + sum( + frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + ) + ) + + # flow_mol_phase_comp Expression line 133 + assert isinstance(frame.props[1].flow_mol_phase_comp, Expression) + assert len(frame.props[1].flow_mol_phase_comp) == 9 + for p, i in frame.props[1].flow_mol_phase_comp: + assert (p, i) in frame.props[1].params._phase_component_set + assert str(frame.props[1].flow_mol_phase_comp[p, i].expr) == str( + frame.props[1].flow_mol_phase[p] + * frame.props[1].mole_frac_phase_comp[p, i] + ) + + @pytest.mark.unit + def test_constraints(self, frame): + # Check that the correct constraints are present + + # mole_frac_comp_eq Constraint line 183 + assert isinstance(frame.props[1].mole_frac_comp_eq, Constraint) + assert len(frame.props[1].mole_frac_comp_eq) == 3 + for i in frame.props[1].mole_frac_comp_eq: + assert str(frame.props[1].mole_frac_comp_eq[i].body) == str( + frame.props[1].mole_frac_comp[i] * frame.props[1].flow_mol + - sum( + frame.props[1].mole_frac_phase_comp[p, i] + * frame.props[1].flow_mol_phase[p] + for p in frame.props[1].params.phase_list + if (p, i) in frame.props[1].params._phase_component_set + ) + ) + + # phase_fraction_constraint Constraint line 200 + assert isinstance(frame.props[1].phase_fraction_constraint, Constraint) + assert len(frame.props[1].phase_fraction_constraint) == 3 + for p in frame.props[1].phase_fraction_constraint: + assert p in frame.params.phase_list + assert str(frame.props[1].phase_fraction_constraint[p].body) == str( + frame.props[1].phase_frac[p] * frame.props[1].flow_mol + - frame.props[1].flow_mol_phase[p] + ) + + assert_units_consistent(frame.props[1]) + + +class TestCommon(object): + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + m.params = DummyParameterBlock( + components={"c1": {}, "c2": {}, "c3": {}}, + phases={ + "a": {"equation_of_state": DummyEoS}, + "b": {"equation_of_state": DummyEoS}, + }, + state_definition=FpTPxpc, + pressure_ref=100000.0, + temperature_ref=300, + state_bounds={ + "flow_mol_phase": (0, 0.1, 0.2, pyunits.kmol / pyunits.s), + "temperature": (522, 621, 720, pyunits.degR), + "pressure": (1, 3, 5, pyunits.bar), + }, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + ) + + # Build state block + m.props = m.params.build_state_block([1], defined_state=True) + + # Add necessary variables that would be built by other methods + m.props[1].dens_mol_phase = Var(m.params.phase_list, initialize=1) + m.props[1].enth_mol_phase = Var(m.params.phase_list, initialize=1) + + return m + + @pytest.mark.unit + def test_convert_vars(self, frame): + # Check that all state var values and bounds were converted correctly + for p in frame.props[1].flow_mol_phase: + assert frame.props[1].flow_mol_phase[p].value == 100 + assert frame.props[1].flow_mol_phase[p].lb == 0 + assert frame.props[1].flow_mol_phase[p].ub == 200 + assert check_units_equivalent( + frame.props[1].flow_mol_phase, pyunits.mol / pyunits.s + ) + + assert frame.props[1].pressure.value == 3e5 + assert frame.props[1].pressure.lb == 1e5 + assert frame.props[1].pressure.ub == 5e5 + assert check_units_equivalent(frame.props[1].pressure, pyunits.Pa) + + assert frame.props[1].temperature.value == 345 + assert frame.props[1].temperature.lb == 290 + assert frame.props[1].temperature.ub == 400 + assert check_units_equivalent(frame.props[1].temperature, pyunits.K) + + # Check supporting variables + + assert isinstance(frame.props[1].flow_mol_phase_comp, Expression) + assert len(frame.props[1].flow_mol_phase_comp) == 6 + + assert isinstance(frame.props[1].mole_frac_comp, Var) + assert len(frame.props[1].mole_frac_comp) == 3 + + assert isinstance(frame.props[1].phase_frac, Var) + assert len(frame.props[1].phase_frac) == 2 + + assert isinstance(frame.props[1].mole_frac_comp_eq, Constraint) + assert len(frame.props[1].mole_frac_comp_eq) == 3 + + assert isinstance(frame.props[1].phase_fraction_constraint, Constraint) + assert len(frame.props[1].phase_fraction_constraint) == 2 + + @pytest.mark.unit + def test_calculate_scaling_factors(self, frame): + frame.props[1].calculate_scaling_factors() + + print(len(frame.props[1].scaling_factor)) + + assert len(frame.props[1].scaling_factor) == 22 + assert frame.props[1].scaling_factor[frame.props[1].flow_mol] == 1e-2 + assert frame.props[1].scaling_factor[frame.props[1].flow_mol_phase["a"]] == 1e-2 + assert frame.props[1].scaling_factor[frame.props[1].flow_mol_phase["b"]] == 1e-2 + assert ( + frame.props[1].scaling_factor[frame.props[1].flow_mol_phase_comp["a", "c1"]] + == 1e-2 + ) + assert ( + frame.props[1].scaling_factor[frame.props[1].flow_mol_phase_comp["a", "c2"]] + == 1e-2 + ) + assert ( + frame.props[1].scaling_factor[frame.props[1].flow_mol_phase_comp["a", "c3"]] + == 1e-2 + ) + assert ( + frame.props[1].scaling_factor[frame.props[1].flow_mol_phase_comp["b", "c1"]] + == 1e-2 + ) + assert ( + frame.props[1].scaling_factor[frame.props[1].flow_mol_phase_comp["b", "c2"]] + == 1e-2 + ) + assert ( + frame.props[1].scaling_factor[frame.props[1].flow_mol_phase_comp["b", "c3"]] + == 1e-2 + ) + assert frame.props[1].scaling_factor[frame.props[1].dens_mol_phase["a"]] == 1e-2 + assert frame.props[1].scaling_factor[frame.props[1].dens_mol_phase["b"]] == 1e-2 + assert frame.props[1].scaling_factor[frame.props[1].mole_frac_comp["c1"]] == 1e3 + assert frame.props[1].scaling_factor[frame.props[1].mole_frac_comp["c2"]] == 1e3 + assert frame.props[1].scaling_factor[frame.props[1].mole_frac_comp["c3"]] == 1e3 + assert ( + frame.props[1].scaling_factor[ + frame.props[1].mole_frac_phase_comp["a", "c1"] + ] + == 1e3 + ) + assert ( + frame.props[1].scaling_factor[ + frame.props[1].mole_frac_phase_comp["a", "c2"] + ] + == 1e3 + ) + assert ( + frame.props[1].scaling_factor[ + frame.props[1].mole_frac_phase_comp["a", "c3"] + ] + == 1e3 + ) + assert ( + frame.props[1].scaling_factor[ + frame.props[1].mole_frac_phase_comp["b", "c1"] + ] + == 1e3 + ) + assert ( + frame.props[1].scaling_factor[ + frame.props[1].mole_frac_phase_comp["b", "c2"] + ] + == 1e3 + ) + assert ( + frame.props[1].scaling_factor[ + frame.props[1].mole_frac_phase_comp["b", "c3"] + ] + == 1e3 + ) + assert frame.props[1].scaling_factor[frame.props[1].pressure] == 1e-5 + assert frame.props[1].scaling_factor[frame.props[1].temperature] == 1e-2 + + # Test General Methods + @pytest.mark.unit + def test_get_material_flow_terms(self, frame): + for p, j in frame.params._phase_component_set: + assert str(frame.props[1].get_material_flow_terms(p, j)) == str( + frame.props[1].flow_mol_phase_comp[p, j] + ) + + @pytest.mark.unit + def test_get_enthalpy_flow_terms(self, frame): + for p in frame.params.phase_list: + assert str(frame.props[1].get_enthalpy_flow_terms(p)) == str( + frame.props[1]._enthalpy_flow_term[p] + ) + assert str(frame.props[1]._enthalpy_flow_term[p].expr) == str( + frame.props[1].flow_mol_phase[p] * frame.props[1].enth_mol_phase[p] + ) + + @pytest.mark.unit + def test_get_material_density_terms(self, frame): + for p, j in frame.params._phase_component_set: + assert str(frame.props[1].get_material_density_terms(p, j)) == str( + frame.props[1]._material_density_term[p, j] + ) + assert str(frame.props[1]._material_density_term[p, j].expr) == str( + frame.props[1].dens_mol_phase[p] + * frame.props[1].mole_frac_phase_comp[p, j] + ) + + @pytest.mark.unit + def test_get_energy_density_terms(self, frame): + for p in frame.params.phase_list: + assert str(frame.props[1].get_energy_density_terms(p)) == str( + frame.props[1]._energy_density_term[p] + ) + assert str(frame.props[1]._energy_density_term[p].expr) == str( + frame.props[1].dens_mol_phase[p] + * frame.props[1].energy_internal_mol_phase[p] + ) + + @pytest.mark.unit + def test_default_material_balance_type(self, frame): + assert ( + frame.props[1].default_material_balance_type() + == MaterialBalanceType.componentTotal + ) + + @pytest.mark.unit + def test_default_energy_balance_type(self, frame): + assert ( + frame.props[1].default_energy_balance_type() + == EnergyBalanceType.enthalpyTotal + ) + + @pytest.mark.unit + def test_get_material_flow_basis(self, frame): + assert frame.props[1].get_material_flow_basis() == MaterialFlowBasis.molar + + @pytest.mark.unit + def test_define_state_vars(self, frame): + assert frame.props[1].define_state_vars() == { + "flow_mol_phase": frame.props[1].flow_mol_phase, + "mole_frac_phase_comp": frame.props[1].mole_frac_phase_comp, + "temperature": frame.props[1].temperature, + "pressure": frame.props[1].pressure, + } + + @pytest.mark.unit + def test_define_display_vars(self, frame): + assert frame.props[1].define_display_vars() == { + "Total Molar Flowrate": frame.props[1].flow_mol_phase, + "Total Mole Fraction": frame.props[1].mole_frac_phase_comp, + "Temperature": frame.props[1].temperature, + "Pressure": frame.props[1].pressure, + } + + @pytest.mark.unit + def test_cloning(self, frame): + """ + This function tests modular properties are cloned correctly using pyomo's + clone method. In particular, it tests that all the methods point to variables + on the clone, and not on the original model. + """ + blk = frame.clone() + + # Test get_material_flow_terms method + for p, j in blk.params._phase_component_set: + assert ( + blk.props[1].get_material_flow_terms(p, j).parent_block() + is blk.props[1] + ) + + # Test get_enthalpy_flow_terms method + for p in blk.params.phase_list: + assert ( + blk.props[1].get_enthalpy_flow_terms(p).parent_block() is blk.props[1] + ) + + # Test get_material_density_terms + for p, j in blk.params._phase_component_set: + assert ( + blk.props[1].get_material_density_terms(p, j).parent_block() + is blk.props[1] + ) + + # Test get_energy_Density_terms + for p in blk.params.phase_list: + assert ( + blk.props[1].get_energy_density_terms(p).parent_block() is blk.props[1] + ) + + for i in blk.props[1].define_state_vars().values(): + assert i.parent_block() is blk.props[1] + + for i in blk.props[1].define_display_vars().values(): + assert i.parent_block() is blk.props[1] + + +# ----------------------------------------------------------------------------- +# Test example from Austin Ladshaw via WaterTAP project +# This example revealed some issues with initialization of FpcTP states with +# phase equilibrium +thermo_config_no_rxn = { + "components": { + "H2O": { + "type": Solvent, + # Define the methods used to calculate the following properties + "dens_mol_liq_comp": Perrys, + "enth_mol_liq_comp": Perrys, + "cp_mol_liq_comp": Perrys, + "entr_mol_liq_comp": Perrys, + "enth_mol_ig_comp": NIST, + "pressure_sat_comp": NIST, + "phase_equilibrium_form": {("Vap", "Liq"): fugacity}, + # Parameter data is always associated with the methods defined above + "parameter_data": { + "mw": (18.0153, pyunits.g / pyunits.mol), + "pressure_crit": (220.64e5, pyunits.Pa), + "temperature_crit": (647, pyunits.K), + # Comes from Perry's Handbook: p. 2-98 + "dens_mol_liq_comp_coeff": { + "eqn_type": 1, + "1": (5.459, pyunits.kmol * pyunits.m**-3), + "2": (0.30542, pyunits.dimensionless), + "3": (647.13, pyunits.K), + "4": (0.081, pyunits.dimensionless), + }, + "enth_mol_form_liq_comp_ref": (-285.830, pyunits.kJ / pyunits.mol), + "enth_mol_form_vap_comp_ref": (0, pyunits.kJ / pyunits.mol), + # Comes from Perry's Handbook: p. 2-174 + "cp_mol_liq_comp_coeff": { + "1": (2.7637e5, pyunits.J / pyunits.kmol / pyunits.K), + "2": (-2.0901e3, pyunits.J / pyunits.kmol / pyunits.K**2), + "3": (8.125, pyunits.J / pyunits.kmol / pyunits.K**3), + "4": (-1.4116e-2, pyunits.J / pyunits.kmol / pyunits.K**4), + "5": (9.3701e-6, pyunits.J / pyunits.kmol / pyunits.K**5), + }, + "cp_mol_ig_comp_coeff": { + "A": (30.09200, pyunits.J / pyunits.mol / pyunits.K), + "B": ( + 6.832514, + pyunits.J * pyunits.mol**-1 * pyunits.K**-1 * pyunits.kiloK**-1, + ), + "C": ( + 6.793435, + pyunits.J * pyunits.mol**-1 * pyunits.K**-1 * pyunits.kiloK**-2, + ), + "D": ( + -2.534480, + pyunits.J * pyunits.mol**-1 * pyunits.K**-1 * pyunits.kiloK**-3, + ), + "E": ( + 0.082139, + pyunits.J * pyunits.mol**-1 * pyunits.K**-1 * pyunits.kiloK**2, + ), + "F": (-250.8810, pyunits.kJ / pyunits.mol), + "G": (223.3967, pyunits.J / pyunits.mol / pyunits.K), + "H": (0, pyunits.kJ / pyunits.mol), + }, + "entr_mol_form_liq_comp_ref": ( + 69.95, + pyunits.J / pyunits.K / pyunits.mol, + ), + "pressure_sat_comp_coeff": { + "A": (4.6543, None), # [1], temperature range 255.9 K - 373 K + "B": (1435.264, pyunits.K), + "C": (-64.848, pyunits.K), + }, + }, + }, + "CO2": { + "type": Solute, + "dens_mol_liq_comp": Perrys, + "enth_mol_liq_comp": Perrys, + "cp_mol_liq_comp": Perrys, + "entr_mol_liq_comp": Perrys, + "enth_mol_ig_comp": NIST, + "pressure_sat_comp": NIST, + "phase_equilibrium_form": {("Vap", "Liq"): fugacity}, + "parameter_data": { + "mw": (44.0095, pyunits.g / pyunits.mol), + "pressure_crit": (73.825e5, pyunits.Pa), + "temperature_crit": (304.23, pyunits.K), + "dens_mol_liq_comp_coeff": { + "eqn_type": 1, + "1": (0.000789, pyunits.kmol * pyunits.m**-3), + "2": (0.000956, pyunits.dimensionless), + "3": (500.78, pyunits.K), + "4": (0.94599, pyunits.dimensionless), + }, + "cp_mol_ig_comp_coeff": { + "A": (24.99735, pyunits.J / pyunits.mol / pyunits.K), + "B": ( + 55.18696, + pyunits.J * pyunits.mol**-1 * pyunits.K**-1 * pyunits.kiloK**-1, + ), + "C": ( + -33.69137, + pyunits.J * pyunits.mol**-1 * pyunits.K**-1 * pyunits.kiloK**-2, + ), + "D": ( + 7.948387, + pyunits.J * pyunits.mol**-1 * pyunits.K**-1 * pyunits.kiloK**-3, + ), + "E": ( + -0.136638, + pyunits.J * pyunits.mol**-1 * pyunits.K**-1 * pyunits.kiloK**2, + ), + "F": (-403.6075, pyunits.kJ / pyunits.mol), + "G": (228.2431, pyunits.J / pyunits.mol / pyunits.K), + "H": (0, pyunits.kJ / pyunits.mol), + }, + "cp_mol_liq_comp_coeff": { + "1": (-8.3043e6, pyunits.J / pyunits.kmol / pyunits.K), + "2": (1.0437e5, pyunits.J / pyunits.kmol / pyunits.K**2), + "3": (4.333e2, pyunits.J / pyunits.kmol / pyunits.K**3), + "4": (6.0052e-1, pyunits.J / pyunits.kmol / pyunits.K**4), + "5": (0, pyunits.J / pyunits.kmol / pyunits.K**5), + }, + "enth_mol_form_liq_comp_ref": (-285.83, pyunits.kJ / pyunits.mol), + "enth_mol_form_vap_comp_ref": (-393.52, pyunits.kJ / pyunits.mol), + "entr_mol_form_liq_comp_ref": (0, pyunits.J / pyunits.K / pyunits.mol), + "entr_mol_form_vap_comp_ref": (213.6, pyunits.J / pyunits.mol), + "pressure_sat_comp_coeff": { + "A": (6.81228, None), + "B": (1301.679, pyunits.K), + "C": (-3.494, pyunits.K), + }, + }, + }, + }, + "phases": { + "Liq": {"type": AqueousPhase, "equation_of_state": Ideal}, + "Vap": {"type": VaporPhase, "equation_of_state": Ideal}, + }, + # Set base units of measurement + "base_units": { + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + # Specifying state definition + "state_definition": FpTPxpc, + "state_bounds": { + "flow_mol_phase": (0, 100, 1000, pyunits.mol / pyunits.s), + "temperature": (273.15, 300, 500, pyunits.K), + "pressure": (5e4, 1e5, 1e6, pyunits.Pa), + }, + "pressure_ref": (101325, pyunits.Pa), + "temperature_ref": (300, pyunits.K), + # Defining phase equilibria + "phases_in_equilibrium": [("Vap", "Liq")], + "phase_equilibrium_state": {("Vap", "Liq"): SmoothVLE}, + "bubble_dew_method": IdealBubbleDew, +} + + +@pytest.mark.component +def test_phase_equilibrium_initializer_object(): + # Create a pyomo model object + model = ConcreteModel() + model.fs = FlowsheetBlock(dynamic=False) + + model.fs.thermo_params = GenericParameterBlock(**thermo_config_no_rxn) + + model.fs.state = model.fs.thermo_params.build_state_block( + model.fs.time, defined_state=False, + ) + + model.fs.state[0].pressure.set_value(101325.0) + model.fs.state[0].temperature.set_value(298.0) + + F_vap = 0.005 + F_liq = 9.995 + + model.fs.state[0].flow_mol_phase["Vap"].set_value(F_vap) + model.fs.state[0].flow_mol_phase["Liq"].set_value(F_liq) + model.fs.state[0].mole_frac_phase_comp["Vap", "CO2"].set_value(0.005 / F_vap) + model.fs.state[0].mole_frac_phase_comp["Vap", "H2O"].set_value(1e-8 / F_vap) + model.fs.state[0].mole_frac_phase_comp["Liq", "CO2"].set_value(1e-8 / F_liq) + model.fs.state[0].mole_frac_phase_comp["Liq", "H2O"].set_value(9.995 / F_liq) + + assert_units_consistent(model) + # We expect 8 state variables, but four additional constraints for phase equilibrium + + assert degrees_of_freedom(model) == 8 - 4 + + initializer = model.fs.state.default_initializer() + + initializer.initialize(model.fs.state) + + # The initialization routine updates the tolerance to handle unconverged constraints + assert initializer.config.constraint_tolerance == float("inf") + + # Check that degrees of freedom are still the same + assert degrees_of_freedom(model) == 8 - 4 + # + # m = model.fs.state[0] + # + # m.display() + # # + # # for v in m.component_data_objects(Var, active=True): + # # print(f"{v.name}: value={v.value}, lb={v.lb}, ub={v.ub}") + # # + # # print() + # # + # # for c in m.component_data_objects(Constraint, active=True): + # # print(c.name, c.expr) + # + # # As the phase equilibrium constraints were not solved, we expect these to have a large residual + # large_res = large_residuals_set(model.fs.state[0]) + # print(large_res) + # assert len(large_res) == 4 + # for i in large_res: + # assert i.name in [ + # "fs.state[0.0].phase_fraction_constraint[Liq]", + # "fs.state[0.0].phase_fraction_constraint[Vap]", + # "fs.state[0.0].equilibrium_constraint[Vap,Liq,H2O]", + # "fs.state[0.0].equilibrium_constraint[Vap,Liq,CO2]", + # ] diff --git a/idaes/models/properties/modular_properties/state_definitions/tests/test_FpTPxpc_electrolyte.py b/idaes/models/properties/modular_properties/state_definitions/tests/test_FpTPxpc_electrolyte.py new file mode 100644 index 0000000000..06a3ea9aa2 --- /dev/null +++ b/idaes/models/properties/modular_properties/state_definitions/tests/test_FpTPxpc_electrolyte.py @@ -0,0 +1,739 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2024 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Tests for constructing and using component lists in electrolyte systems for FpTPxpc state formulation + +Authors: Tanner Polley +""" +# Import Python libraries +from copy import deepcopy +import pytest + +# Import Pyomo units +from pyomo.environ import ( + check_optimal_termination, + ConcreteModel, + Constraint, + Expression, + Set, + value, + Var, + units as pyunits, +) + +# Import IDAES cores +from idaes.models.properties.modular_properties.state_definitions import FpTPxpc +from idaes.models.properties.modular_properties.state_definitions.tests.electrolyte_testing_config import ( + get_config_no_inherent_reactions, + get_config_with_inherent_reactions, +) +from idaes.core import FlowsheetBlock +from idaes.models.properties.modular_properties.base.generic_property import ( + GenericParameterBlock, + StateIndex, +) + +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.solvers import get_solver + +solver = get_solver("ipopt_v2") + +state_bounds = { + "flow_mol_phase": (0, 100, 1000, pyunits.mol / pyunits.s), + "temperature": (273.15, 300, 500, pyunits.K), + "pressure": (5e4, 1e5, 1e6, pyunits.Pa), +} + + +# ----------------------------------------------------------------------------- +class TestApparentSpeciesBasisNoInherent: + config = get_config_no_inherent_reactions() + config["state_components"] = StateIndex.apparent + config["state_definition"] = FpTPxpc + config["state_bounds"] = deepcopy(state_bounds) + + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = GenericParameterBlock(**TestApparentSpeciesBasisNoInherent.config) + + m.fs.state = m.fs.props.build_state_block([1], defined_state=True) + + return m + + @pytest.mark.unit + def test_vars_and_constraints(self, frame): + m = frame + + assert m.fs.state[1].component_list is m.fs.props.apparent_species_set + + assert ( + m.fs.state[1].phase_component_set is m.fs.props.apparent_phase_component_set + ) + + assert isinstance(m.fs.state[1].flow_mol_phase_comp, Expression) + assert len(m.fs.state[1].flow_mol_phase_comp) == 7 + assert isinstance(m.fs.state[1].pressure, Var) + assert len(m.fs.state[1].pressure) == 1 + assert isinstance(m.fs.state[1].temperature, Var) + assert len(m.fs.state[1].temperature) == 1 + + assert isinstance(m.fs.state[1].flow_mol_phase, Var) + assert len(m.fs.state[1].flow_mol_phase) == 2 + for p in m.fs.state[1].flow_mol_phase: + assert p in ["Liq", "Vap"] + assert isinstance(m.fs.state[1].phase_frac, Var) + assert len(m.fs.state[1].phase_frac) == 2 + for p in m.fs.state[1].phase_frac: + assert p in ["Liq", "Vap"] + + assert isinstance(m.fs.state[1].mole_frac_comp, Var) + assert len(m.fs.state[1].mole_frac_comp) == 4 + for j in m.fs.state[1].mole_frac_comp: + assert j in ["KHCO3", "H2O", "CO2", "N2"] + + assert isinstance(m.fs.state[1].mole_frac_phase_comp, Var) + assert len(m.fs.state[1].mole_frac_phase_comp) == 7 + for j in m.fs.state[1].mole_frac_phase_comp: + assert j in [ + ("Liq", "H2O"), + ("Liq", "CO2"), + ("Liq", "KHCO3"), + ("Vap", "H2O"), + ("Vap", "CO2"), + ("Vap", "KHCO3"), + ("Vap", "N2"), + ] + + # Check references to base state variables + assert m.fs.state[1].flow_mol_apparent is m.fs.state[1].flow_mol + for i in m.fs.state[1].flow_mol_phase_apparent: + assert ( + m.fs.state[1].flow_mol_phase_apparent[i] + is m.fs.state[1].flow_mol_phase[i] + ) + assert i in m.fs.props.phase_list + for i in m.fs.state[1].flow_mol_phase_comp_apparent: + assert ( + m.fs.state[1].flow_mol_phase_comp_apparent[i] + is m.fs.state[1].flow_mol_phase_comp[i] + ) + assert i in m.fs.props.apparent_phase_component_set + for i in m.fs.state[1].mole_frac_phase_comp_apparent: + assert ( + m.fs.state[1].mole_frac_phase_comp_apparent[i] + is m.fs.state[1].mole_frac_phase_comp[i] + ) + assert i in m.fs.props.apparent_phase_component_set + + # Check for true species components + assert isinstance(m.fs.state[1].flow_mol_phase_comp_true, Var) + assert len(m.fs.state[1].flow_mol_phase_comp_true) == 8 + assert isinstance(m.fs.state[1].appr_to_true_species, Constraint) + assert len(m.fs.state[1].appr_to_true_species) == 8 + + @pytest.mark.component + def test_solve_for_true_species(self, frame): + m = frame + + m.fs.state[1].flow_mol_phase["Liq"].fix(1.0) + m.fs.state[1].flow_mol_phase["Vap"].fix(0.5 + (1 / 6) * 3) + m.fs.state[1].mole_frac_phase_comp["Liq", "H2O"].fix( + 1 / 3 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "CO2"].fix( + 1 / 3 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "KHCO3"].fix( + 1 / 3 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Vap", "H2O"].fix( + 1 / 6 / value(m.fs.state[1].flow_mol_phase["Vap"]) + ) + m.fs.state[1].mole_frac_phase_comp["Vap", "CO2"].fix( + 1 / 6 / value(m.fs.state[1].flow_mol_phase["Vap"]) + ) + m.fs.state[1].mole_frac_phase_comp["Vap", "KHCO3"].fix( + 1 / 6 / value(m.fs.state[1].flow_mol_phase["Vap"]) + ) + m.fs.state[1].mole_frac_phase_comp["Vap", "N2"].fix( + 0.5 / value(m.fs.state[1].flow_mol_phase["Vap"]) + ) + m.fs.state[1].temperature.fix(300) + m.fs.state[1].pressure.fix(1e5) + + assert degrees_of_freedom(m.fs) == 0 + + res = solver.solve(m.fs, tee=True) + + # Check for optimal solution + assert check_optimal_termination(res) + + # Check true species flowrates + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Vap", "H2O"] + ) == pytest.approx(1 / 6, rel=1e-5) + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Vap", "CO2"] + ) == pytest.approx(1 / 6, rel=1e-5) + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Vap", "KHCO3"] + ) == pytest.approx(1 / 6, rel=1e-5) + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Vap", "N2"] + ) == pytest.approx(0.5, rel=1e-5) + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "H2O"] + ) == pytest.approx(1 / 3, rel=1e-5) + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "CO2"] + ) == pytest.approx(1 / 3, rel=1e-5) + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "K+"] + ) == pytest.approx(1 / 3, rel=1e-5) + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "HCO3-"] + ) == pytest.approx(1 / 3, rel=1e-5) + + +# ----------------------------------------------------------------------------- +def dens_mol_H2O(*args, **kwargs): + return 55e3 + + +class TestApparentSpeciesBasisInherent: + config = get_config_with_inherent_reactions() + config["state_components"] = StateIndex.apparent + config["state_definition"] = FpTPxpc + config["state_bounds"] = deepcopy(state_bounds) + + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = GenericParameterBlock(**TestApparentSpeciesBasisInherent.config) + + m.fs.state = m.fs.props.build_state_block([1], defined_state=True) + + m.fs.state[1].calculate_scaling_factors() + + return m + + @pytest.mark.unit + def test_vars_and_constraints(self, frame): + m = frame + + assert m.fs.state[1].component_list is m.fs.props.apparent_species_set + + assert ( + m.fs.state[1].phase_component_set is m.fs.props.apparent_phase_component_set + ) + + assert isinstance(m.fs.state[1].flow_mol_phase_comp, Expression) + assert len(m.fs.state[1].flow_mol_phase_comp) == 4 + assert isinstance(m.fs.state[1].pressure, Var) + assert len(m.fs.state[1].pressure) == 1 + assert isinstance(m.fs.state[1].temperature, Var) + assert len(m.fs.state[1].temperature) == 1 + + assert isinstance(m.fs.state[1].flow_mol_phase, Var) + assert len(m.fs.state[1].flow_mol_phase) == 1 + for p in m.fs.state[1].flow_mol_phase: + assert p in ["Liq"] + assert isinstance(m.fs.state[1].phase_frac, Var) + assert len(m.fs.state[1].phase_frac) == 1 + for p in m.fs.state[1].phase_frac: + assert p in ["Liq"] + + assert isinstance(m.fs.state[1].mole_frac_comp, Var) + assert len(m.fs.state[1].mole_frac_comp) == 4 + for j in m.fs.state[1].mole_frac_comp: + assert j in ["KHCO3", "H2O", "K2CO3", "KOH"] + + assert isinstance(m.fs.state[1].mole_frac_phase_comp, Var) + assert len(m.fs.state[1].mole_frac_phase_comp) == 4 + for j in m.fs.state[1].mole_frac_phase_comp: + assert j in [ + ("Liq", "H2O"), + ("Liq", "K2CO3"), + ("Liq", "KHCO3"), + ("Liq", "KOH"), + ] + + # Check references to base state variables + assert m.fs.state[1].flow_mol_apparent is m.fs.state[1].flow_mol + for i in m.fs.state[1].flow_mol_phase_apparent: + assert ( + m.fs.state[1].flow_mol_phase_apparent[i] + is m.fs.state[1].flow_mol_phase[i] + ) + assert i in m.fs.props.phase_list + for i in m.fs.state[1].flow_mol_phase_comp_apparent: + assert ( + m.fs.state[1].flow_mol_phase_comp_apparent[i] + is m.fs.state[1].flow_mol_phase_comp[i] + ) + assert i in m.fs.props.apparent_phase_component_set + for i in m.fs.state[1].mole_frac_phase_comp_apparent: + assert ( + m.fs.state[1].mole_frac_phase_comp_apparent[i] + is m.fs.state[1].mole_frac_phase_comp[i] + ) + assert i in m.fs.props.apparent_phase_component_set + + # Check for true species components + assert isinstance(m.fs.state[1].flow_mol_phase_comp_true, Var) + assert len(m.fs.state[1].flow_mol_phase_comp_true) == 6 + assert isinstance(m.fs.state[1].mole_frac_phase_comp_true, Var) + assert len(m.fs.state[1].mole_frac_phase_comp_true) == 6 + assert isinstance(m.fs.state[1].appr_to_true_species, Constraint) + assert len(m.fs.state[1].appr_to_true_species) == 6 + assert isinstance(m.fs.state[1].true_mole_frac_constraint, Constraint) + assert len(m.fs.state[1].true_mole_frac_constraint) == 6 + + # Check for inherent reactions + assert m.fs.state[1].has_inherent_reactions + assert not m.fs.state[1].include_inherent_reactions + assert isinstance(m.fs.state[1].params.inherent_reaction_idx, Set) + assert len(m.fs.state[1].params.inherent_reaction_idx) == 2 + for i in m.fs.state[1].params.inherent_reaction_idx: + assert i in ["h2o_si", "co3_hco3"] + + assert isinstance(m.fs.state[1].apparent_inherent_reaction_extent, Var) + assert len(m.fs.state[1].apparent_inherent_reaction_extent) == 2 + + @pytest.mark.component + def test_solve_for_true_species(self, frame): + m = frame + + m.fs.state[1].flow_mol_phase["Liq"].fix(1.6 + 0.4 + 1e-8 + 1e-8) + m.fs.state[1].mole_frac_phase_comp["Liq", "H2O"].fix( + 1.6 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "K2CO3"].fix( + 0.4 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "KHCO3"].fix( + 1e-8 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "KOH"].fix( + 1e-8 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].temperature.fix(350) + m.fs.state[1].pressure.fix(1e5) + + assert degrees_of_freedom(m.fs) == 0 + + m.fs.state.initialize() + + res = solver.solve(m.fs) + + # Check for optimal solution + assert check_optimal_termination(res) + + # Check apparent species flowrates + for j in m.fs.state[1].mole_frac_comp: + print() + print(j) + print(value(m.fs.state[1].flow_mol_phase_comp_apparent["Liq", j])) + print(value(m.fs.state[1].flow_mol), value(m.fs.state[1].mole_frac_comp[j])) + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", j] + ) == pytest.approx( + value(m.fs.state[1].flow_mol * m.fs.state[1].mole_frac_comp[j]), + rel=1e-5, + ) + + # Check element balances + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "K2CO3"] * 2 + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KHCO3"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KOH"] + ) == pytest.approx( + value(m.fs.state[1].flow_mol_phase_comp_true["Liq", "K+"]), rel=1e-5 + ) + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "K2CO3"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KHCO3"] + ) == pytest.approx( + value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "CO3--"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "HCO3-"] + ), + rel=1e-5, + ) + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "K2CO3"] * 3 + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KHCO3"] * 3 + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KOH"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "H2O"] + ) == pytest.approx( + value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "CO3--"] * 3 + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "HCO3-"] * 3 + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "OH-"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "H2O"] + ), + rel=1e-5, + ) + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KHCO3"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KOH"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "H2O"] * 2 + ) == pytest.approx( + value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "HCO3-"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "OH-"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "H+"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "H2O"] * 2 + ), + rel=1e-5, + ) + + # Check true species mole fractions + assert value( + m.fs.state[1].mole_frac_phase_comp_true["Liq", "CO3--"] + ) == pytest.approx(0.142857, rel=1e-5) + assert value( + m.fs.state[1].mole_frac_phase_comp_true["Liq", "H+"] + ) == pytest.approx(2.90081e-16, rel=1e-5) + assert value( + m.fs.state[1].mole_frac_phase_comp_true["Liq", "H2O"] + ) == pytest.approx(0.571429, rel=1e-5) + assert value( + m.fs.state[1].mole_frac_phase_comp_true["Liq", "HCO3-"] + ) == pytest.approx(1.13961e-08, rel=1e-5) + assert value( + m.fs.state[1].mole_frac_phase_comp_true["Liq", "K+"] + ) == pytest.approx(0.285714, rel=1e-5) + assert value( + m.fs.state[1].mole_frac_phase_comp_true["Liq", "OH-"] + ) == pytest.approx(1.139606e-08, rel=1e-5) + + +# ----------------------------------------------------------------------------- +class TestTrueSpeciesBasisNoInherent: + config = get_config_no_inherent_reactions() + config["state_components"] = StateIndex.true + config["state_definition"] = FpTPxpc + config["state_bounds"] = deepcopy(state_bounds) + + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = GenericParameterBlock(**TestTrueSpeciesBasisNoInherent.config) + + m.fs.state = m.fs.props.build_state_block([1], defined_state=True) + + return m + + @pytest.mark.unit + def test_vars_and_constraints(self, frame): + m = frame + + assert m.fs.state[1].component_list is m.fs.props.true_species_set + + assert m.fs.state[1].phase_component_set is m.fs.props.true_phase_component_set + + assert isinstance(m.fs.state[1].flow_mol_phase_comp, Expression) + assert len(m.fs.state[1].flow_mol_phase_comp) == 8 + assert isinstance(m.fs.state[1].pressure, Var) + assert len(m.fs.state[1].pressure) == 1 + assert isinstance(m.fs.state[1].temperature, Var) + assert len(m.fs.state[1].temperature) == 1 + + assert isinstance(m.fs.state[1].flow_mol_phase, Var) + assert len(m.fs.state[1].flow_mol_phase) == 2 + for p in m.fs.state[1].flow_mol_phase: + assert p in ["Liq", "Vap"] + assert isinstance(m.fs.state[1].phase_frac, Var) + assert len(m.fs.state[1].phase_frac) == 2 + for p in m.fs.state[1].phase_frac: + assert p in ["Liq", "Vap"] + + assert isinstance(m.fs.state[1].mole_frac_comp, Var) + assert len(m.fs.state[1].mole_frac_comp) == 5 + for j in m.fs.state[1].mole_frac_comp: + assert j in ["K+", "HCO3-", "H2O", "CO2", "N2"] + + assert isinstance(m.fs.state[1].mole_frac_phase_comp, Var) + assert len(m.fs.state[1].mole_frac_phase_comp) == 8 + for j in m.fs.state[1].mole_frac_phase_comp: + assert j in [ + ("Liq", "H2O"), + ("Liq", "CO2"), + ("Liq", "K+"), + ("Liq", "HCO3-"), + ("Vap", "H2O"), + ("Vap", "CO2"), + ("Vap", "KHCO3"), + ("Vap", "N2"), + ] + + # Check references to base state variables + assert m.fs.state[1].flow_mol_true is m.fs.state[1].flow_mol + for i in m.fs.state[1].flow_mol_phase_true: + assert ( + m.fs.state[1].flow_mol_phase_true[i] is m.fs.state[1].flow_mol_phase[i] + ) + assert i in m.fs.props.phase_list + for i in m.fs.state[1].flow_mol_phase_comp_true: + assert ( + m.fs.state[1].flow_mol_phase_comp_true[i] + is m.fs.state[1].flow_mol_phase_comp[i] + ) + assert i in m.fs.props.true_phase_component_set + for i in m.fs.state[1].mole_frac_phase_comp_true: + assert ( + m.fs.state[1].mole_frac_phase_comp_true[i] + is m.fs.state[1].mole_frac_phase_comp[i] + ) + assert i in m.fs.props.true_phase_component_set + + # Check for apparent species components + assert isinstance(m.fs.state[1].flow_mol_phase_comp_apparent, Var) + assert len(m.fs.state[1].flow_mol_phase_comp_apparent) == 7 + assert isinstance(m.fs.state[1].true_to_appr_species, Constraint) + assert len(m.fs.state[1].true_to_appr_species) == 7 + + +# ----------------------------------------------------------------------------- +class TestTrueSpeciesBasisInherent: + config = get_config_with_inherent_reactions() + config["state_components"] = StateIndex.true + config["state_definition"] = FpTPxpc + config["state_bounds"] = deepcopy(state_bounds) + + @pytest.fixture(scope="class") + def frame(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = GenericParameterBlock(**TestTrueSpeciesBasisInherent.config) + + m.fs.state = m.fs.props.build_state_block([1], defined_state=True) + + m.fs.state[1].calculate_scaling_factors() + + return m + + @pytest.mark.unit + def test_vars_and_constraints(self, frame): + m = frame + + assert m.fs.state[1].component_list is m.fs.props.true_species_set + + assert m.fs.state[1].phase_component_set is m.fs.props.true_phase_component_set + + assert isinstance(m.fs.state[1].flow_mol_phase_comp, Expression) + assert len(m.fs.state[1].flow_mol_phase_comp) == 6 + assert isinstance(m.fs.state[1].pressure, Var) + assert len(m.fs.state[1].pressure) == 1 + assert isinstance(m.fs.state[1].temperature, Var) + assert len(m.fs.state[1].temperature) == 1 + + assert isinstance(m.fs.state[1].flow_mol_phase, Var) + assert len(m.fs.state[1].flow_mol_phase) == 1 + for p in m.fs.state[1].flow_mol_phase: + assert p in ["Liq"] + assert isinstance(m.fs.state[1].phase_frac, Var) + assert len(m.fs.state[1].phase_frac) == 1 + for p in m.fs.state[1].phase_frac: + assert p in ["Liq"] + + assert isinstance(m.fs.state[1].mole_frac_comp, Var) + assert len(m.fs.state[1].mole_frac_comp) == 6 + for j in m.fs.state[1].mole_frac_comp: + assert j in ["H+", "K+", "HCO3-", "CO3--", "OH-", "H2O"] + + assert isinstance(m.fs.state[1].mole_frac_phase_comp, Var) + assert len(m.fs.state[1].mole_frac_phase_comp) == 6 + for j in m.fs.state[1].mole_frac_phase_comp: + assert j in [ + ("Liq", "H+"), + ("Liq", "K+"), + ("Liq", "HCO3-"), + ("Liq", "CO3--"), + ("Liq", "OH-"), + ("Liq", "H2O"), + ] + + # Check references to base state variables + assert m.fs.state[1].flow_mol_true is m.fs.state[1].flow_mol + for i in m.fs.state[1].flow_mol_phase_true: + assert ( + m.fs.state[1].flow_mol_phase_true[i] is m.fs.state[1].flow_mol_phase[i] + ) + assert i in m.fs.props.phase_list + for i in m.fs.state[1].flow_mol_phase_comp_true: + assert ( + m.fs.state[1].flow_mol_phase_comp_true[i] + is m.fs.state[1].flow_mol_phase_comp[i] + ) + assert i in m.fs.props.true_phase_component_set + for i in m.fs.state[1].mole_frac_phase_comp_true: + assert ( + m.fs.state[1].mole_frac_phase_comp_true[i] + is m.fs.state[1].mole_frac_phase_comp[i] + ) + assert i in m.fs.props.true_phase_component_set + + # Check for apparent species components + assert isinstance(m.fs.state[1].flow_mol_phase_comp_apparent, Var) + assert len(m.fs.state[1].flow_mol_phase_comp_apparent) == 4 + assert isinstance(m.fs.state[1].mole_frac_phase_comp_apparent, Var) + assert len(m.fs.state[1].mole_frac_phase_comp_apparent) == 4 + assert isinstance(m.fs.state[1].true_to_appr_species, Constraint) + assert len(m.fs.state[1].true_to_appr_species) == 4 + assert isinstance(m.fs.state[1].appr_mole_frac_constraint, Constraint) + assert len(m.fs.state[1].appr_mole_frac_constraint) == 4 + + # Check for inherent reactions + assert m.fs.state[1].has_inherent_reactions + assert m.fs.state[1].include_inherent_reactions + assert isinstance(m.fs.state[1].params.inherent_reaction_idx, Set) + assert len(m.fs.state[1].params.inherent_reaction_idx) == 2 + for i in m.fs.state[1].params.inherent_reaction_idx: + assert i in ["h2o_si", "co3_hco3"] + + assert not hasattr(m.fs.state[1], "apparent_inherent_reaction_extent") + + @pytest.mark.component + def test_solve_for_true_species(self, frame): + m = frame + + m.fs.state[1].temperature.fix(350) + m.fs.state[1].pressure.fix(1e5) + + flow_mol_phase_comp_H2O = 0.5716 * 2 + flow_mol_phase_comp_K = 0.2858 * 2 + flow_mol_phase_comp_H = 2.9e-16 * 2 + flow_mol_phase_comp_HCO3 = 1.14e-8 * 2 + flow_mol_phase_comp_CO3 = 0.1429 * 2 + flow_mol_phase_comp_OH = 1.14e-8 * 2 + + m.fs.state[1].flow_mol_phase["Liq"].fix( + flow_mol_phase_comp_H2O + + flow_mol_phase_comp_K + + flow_mol_phase_comp_H + + flow_mol_phase_comp_HCO3 + + flow_mol_phase_comp_CO3 + + flow_mol_phase_comp_OH + ) + + m.fs.state[1].mole_frac_phase_comp["Liq", "H2O"].fix( + flow_mol_phase_comp_H2O / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "K+"].fix( + flow_mol_phase_comp_K / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "H+"].fix( + flow_mol_phase_comp_H / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "HCO3-"].fix( + flow_mol_phase_comp_HCO3 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "CO3--"].fix( + flow_mol_phase_comp_CO3 / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + m.fs.state[1].mole_frac_phase_comp["Liq", "OH-"].fix( + flow_mol_phase_comp_OH / value(m.fs.state[1].flow_mol_phase["Liq"]) + ) + + assert degrees_of_freedom(m.fs) == 0 + + m.fs.state.initialize() + + res = solver.solve(m.fs) + + # Check for optimal solution + assert check_optimal_termination(res) + + # Check true species flowrates + for j in m.fs.state[1].mole_frac_comp: + assert value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", j] + ) == pytest.approx( + value(m.fs.state[1].flow_mol * m.fs.state[1].mole_frac_comp[j]), + rel=1e-5, + ) + + # Check element balances + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "K2CO3"] * 2 + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KHCO3"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KOH"] + ) == pytest.approx( + value(m.fs.state[1].flow_mol_phase_comp_true["Liq", "K+"]), rel=1e-5 + ) + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "K2CO3"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KHCO3"] + ) == pytest.approx( + value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "CO3--"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "HCO3-"] + ), + rel=1e-5, + ) + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "K2CO3"] * 3 + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KHCO3"] * 3 + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KOH"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "H2O"] + ) == pytest.approx( + value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "CO3--"] * 3 + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "HCO3-"] * 3 + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "OH-"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "H2O"] + ), + rel=1e-5, + ) + assert value( + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KHCO3"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "KOH"] + + m.fs.state[1].flow_mol_phase_comp_apparent["Liq", "H2O"] * 2 + ) == pytest.approx( + value( + m.fs.state[1].flow_mol_phase_comp_true["Liq", "HCO3-"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "OH-"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "H+"] + + m.fs.state[1].flow_mol_phase_comp_true["Liq", "H2O"] * 2 + ), + rel=1e-5, + ) + + # Check apparent species mole fractions + m.fs.state[1].mole_frac_phase_comp_apparent.display() + assert value( + m.fs.state[1].mole_frac_phase_comp_apparent["Liq", "K2CO3"] + ) == pytest.approx(0.2, rel=1e-5) + assert value( + m.fs.state[1].mole_frac_phase_comp_apparent["Liq", "H2O"] + ) == pytest.approx(0.8, rel=1e-5) + assert value( + m.fs.state[1].mole_frac_phase_comp_apparent["Liq", "KHCO3"] + ) == pytest.approx(0, abs=1e-6) + assert value( + m.fs.state[1].mole_frac_phase_comp_apparent["Liq", "KHCO3"] + ) == pytest.approx(0, abs=1e-6)