diff --git a/property_packages/tests/__init__.py b/property_packages/base/__init__.py similarity index 100% rename from property_packages/tests/__init__.py rename to property_packages/base/__init__.py diff --git a/property_packages/base/state_block_constraints.py b/property_packages/base/state_block_constraints.py new file mode 100644 index 0000000..fe0d567 --- /dev/null +++ b/property_packages/base/state_block_constraints.py @@ -0,0 +1,44 @@ +from pyomo.environ import Block, Constraint +from pyomo.core.base.expression import ScalarExpression, Expression, _GeneralExpressionData, ExpressionData +from pyomo.core.base.var import ScalarVar, _GeneralVarData, VarData, IndexedVar, Var + +from property_packages.utils.add_extra_expressions import add_extra_expressions + + +class StateBlockConstraints: + + """ + Extended state blocks should inherit from this class. + Adds additional expressions and constraints to a state block. + """ + + def build(blk, *args): + add_extra_expressions(blk) + blk.constraints = Block() + + + def constrain(blk, name: str, value: float) -> Constraint | Var | None: + """constrain a component by name to a value""" + var = getattr(blk, name) + return blk.constrain_component(var, value) + + + def constrain_component(blk, component: Var | Expression, value: float) -> Constraint | Var | None: + """ + Constrain a component to a value + """ + if type(component) == ScalarExpression: + c = Constraint(expr=component == value) + c.defining_state_var = True + blk.constraints.add_component(component.local_name, c) + return c + elif type(component) in (ScalarVar, _GeneralVarData, VarData, IndexedVar): + component.fix(value) + return component + elif type(component) in (_GeneralExpressionData, ExpressionData): + # allowed, but we don't need to fix it (eg. mole_frac_comp in helmholtz) + return None + else: + raise Exception( + f"Component {component} is not a Var or Expression: {type(component)}" + ) diff --git a/property_packages/helmholtz/helmholtz_extended.py b/property_packages/helmholtz/helmholtz_extended.py index 290f504..d0d6377 100644 --- a/property_packages/helmholtz/helmholtz_extended.py +++ b/property_packages/helmholtz/helmholtz_extended.py @@ -1,101 +1,87 @@ -# extends the IDAES helmholtz property package to include additional properties and methods. +from pyomo.environ import ( + Expression, + SolverFactory, + check_optimal_termination, + units +) from idaes.models.properties.general_helmholtz.helmholtz_state import HelmholtzStateBlockData, _StateBlock from idaes.models.properties.general_helmholtz.helmholtz_functions import HelmholtzParameterBlockData from idaes.core import declare_process_block_class -from property_packages.utils.add_extra_expressions import add_extra_expressions -from pyomo.environ import Constraint, Block -from pyomo.core.base.expression import Expression, ScalarExpression, _GeneralExpressionData, ExpressionData -from pyomo.core.base.var import IndexedVar, ScalarVar, Var, _GeneralVarData,VarData +from idaes.core.util.initialization import solve_indexed_blocks, revert_state_vars +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.exceptions import InitializationError import idaes.logger as idaeslog +from property_packages.base.state_block_constraints import StateBlockConstraints +from property_packages.utils.fix_state_vars import fix_state_vars + class _ExtendedStateBlock(_StateBlock): - """ - This class contains methods which should be applied to Property Blocks as a - whole, rather than individual elements of indexed Property Blocks. - """ - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - def initialize(self, *args, **kwargs): - hold_state = kwargs.pop("hold_state", False) - for i, v in self.items(): - v.constraints.deactivate() - res = super().initialize(*args, **kwargs) - flags = {} - for i, v in self.items(): - v.constraints.activate() - flags[i] = {} - if hold_state: - # Fix the required state variables for zero degrees of freedom, and return a dictionary of the flags. - if not hasattr(v.constraints,"flow_mass") and not v.flow_mol.is_fixed(): - # We need to fix the flow_mol variable - flags[i]["flow_mol"] = True - v.flow_mol.fix() - avaliable_constraints = ["enth_mass","temperature","total_energy_flow","entr_mass","entr_mol","smooth_temperature","custom_vapor_frac"] - if not v.enth_mol.is_fixed(): - # check if any of the constraints exist - found_constraint = False - for constraint in avaliable_constraints: - if hasattr(v.constraints,constraint): - # we don't need to fix the variable - # but we need to remove this from the list of constraints (it can't be used to define pressure) - avaliable_constraints.remove(constraint) - found_constraint = True - break - if not found_constraint: - # we need to fix the variable - flags[i]["enth_mol"] = True - v.enth_mol.fix() - if not v.pressure.is_fixed(): - # check if any of the constraints exist - found_constraint = False - for constraint in avaliable_constraints: - if hasattr(v.constraints,constraint): - # we don't need to fix the variable - # but we need to remove this from the list of constraints (it can't be used to define pressure) - avaliable_constraints.remove(constraint) - found_constraint = True - break - if not found_constraint: - # we need to fix the variable - flags[i]["pressure"] = True - v.pressure.fix() - return flags - - def release_state(self, flags, outlvl=idaeslog.NOTSET): - for i, v in self.items(): - for key in flags[i]: - getattr(v,key).unfix() + def initialize(blk, *args, **kwargs): + outlvl = kwargs.get("outlvl", idaeslog.NOTSET) + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties") + + flag_dict = fix_state_vars(blk, kwargs.get("state_args", None)) + + dof = degrees_of_freedom(blk) + if dof != 0: + raise InitializationError( + f"{blk.name} Unexpected degrees of freedom during " + f"initialization at property initialization step: {dof}." + ) + + res = None + opt = SolverFactory('ipopt') + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + try: + res = solve_indexed_blocks(opt, [blk], tee=slc.tee) + except ValueError as e: + if str(e).startswith("No variables appear"): + # https://github.com/Pyomo/pyomo/pull/3445 + pass + else: + raise e + + if res is not None and not check_optimal_termination(res): + raise InitializationError( + f"{blk.name} failed to initialize successfully. Please check " + f"the output logs for more information." + ) + + if kwargs.get("hold_state") is True: + return flag_dict + else: + blk.release_state(flag_dict) + + init_log.info( + "Property package initialization: {}.".format(idaeslog.condition(res)) + ) + def release_state(blk, flags, outlvl=idaeslog.NOTSET): + revert_state_vars(blk, flags) + @declare_process_block_class("HelmholtzExtendedStateBlock", block_class=_ExtendedStateBlock) -class HelmholtzExtendedStateBlockData(HelmholtzStateBlockData): +class HelmholtzExtendedStateBlockData(HelmholtzStateBlockData, StateBlockConstraints): + + def build(blk, *args): + HelmholtzStateBlockData.build(blk, *args) + StateBlockConstraints.build(blk, *args) + + # rename temperature to old_temperature + old_temperature = blk.temperature + blk.del_component("temperature") + blk.add_component("old_temperature", old_temperature) + # create smooth temperature expression + blk.add_component("temperature", Expression( + expr=blk.old_temperature + (blk.enth_mol / (1 * units.J/units.mol))*0.000001 * units.K + )) - def build(self, *args): - super().build(*args) - # Add expressions for smooth_temperature, enthalpy in terms of mass, etc. - add_extra_expressions(self) - # Add a block for constraints, so we can disable or enable them in bulk - self.constraints = Block() - - def constrain(self,name:str,value:float): - # Value must be a float. TODO: Handle unit conversion. - var = getattr(self,name) - if type(var) == ScalarExpression: - self.constraints.add_component(name, Constraint(expr=var == value)) - elif type(var) in (ScalarVar, _GeneralVarData, VarData): - var.fix(value) - elif type(var) in ( _GeneralExpressionData, ExpressionData) : - # allowed, but we don't need to fix it (eg. mole_frac_comp in helmholtz) - print(f"Variable {self} {name} is an Expression: {type(var)}") - pass - else: - raise Exception(f"Variable {self} {name} is not a Var or Expression: {type(var)}") @declare_process_block_class("HelmholtzExtendedParameterBlock") class HelmholtzExtendedParameterBlockData(HelmholtzParameterBlockData): def build(self): super().build() - # set that we should use the extended state block - self._state_block_class = HelmholtzExtendedStateBlock # type: ignore because it'll get created with declare_process_block_class + self._state_block_class = HelmholtzExtendedStateBlock # noqa: F821 diff --git a/property_packages/tests/test_helmholtz.py b/property_packages/helmholtz/tests/test_helmholtz.py similarity index 83% rename from property_packages/tests/test_helmholtz.py rename to property_packages/helmholtz/tests/test_helmholtz.py index 7ab5359..72359f9 100644 --- a/property_packages/tests/test_helmholtz.py +++ b/property_packages/helmholtz/tests/test_helmholtz.py @@ -1,5 +1,5 @@ # Build and solve a heater block. -from ..build_package import build_package +from property_packages.build_package import build_package from pytest import approx # Import objects from pyomo package @@ -9,8 +9,6 @@ from idaes.core import FlowsheetBlock from idaes.core.util.model_statistics import degrees_of_freedom from idaes.models.unit_models.heater import Heater -from idaes.core.util.tables import _get_state_from_port - def test_helmholtz(): @@ -23,12 +21,12 @@ def test_helmholtz(): m.fs.heater.inlet.flow_mol.fix(1) #m.fs.heater.inlet.vapor_frac.fix(0) #m.fs.heater.inlet.temperature.fix(298) # room temperature in K - m.fs.heater.inlet.enth_mol.fix(1878.87) + m.fs.heater.inlet.enth_mol.fix(1878.71) m.fs.heater.inlet.pressure.fix(101325) assert degrees_of_freedom(m) == 0 m.fs.heater.initialize() solver = SolverFactory('ipopt') solver.solve(m, tee=True) - assert value(_get_state_from_port(m.fs.heater.outlet,0).temperature) == approx(298) + assert value(m.fs.heater.control_volume.properties_out[0].temperature) == approx(298) assert value(m.fs.heater.outlet.pressure[0]) == approx(101325) assert value(m.fs.heater.outlet.flow_mol[0]) == approx(1) \ No newline at end of file diff --git a/property_packages/helmholtz/tests/test_state_definitions.py b/property_packages/helmholtz/tests/test_state_definitions.py index bb58257..4212f25 100644 --- a/property_packages/helmholtz/tests/test_state_definitions.py +++ b/property_packages/helmholtz/tests/test_state_definitions.py @@ -1,8 +1,10 @@ from ..helmholtz_builder import build_helmholtz_package from pytest import approx from pyomo.environ import ConcreteModel, SolverFactory, value, units +from pyomo.core.base.constraint import Constraint, ScalarConstraint from idaes.core import FlowsheetBlock from idaes.models.unit_models import Compressor +from idaes.core.util.model_statistics import degrees_of_freedom def flowsheet(): m = ConcreteModel() @@ -10,30 +12,57 @@ def flowsheet(): m.fs.properties = build_helmholtz_package(["h2o"]) return m + +def build_state(m): + m.fs.state = m.fs.properties.build_state_block([0], defined_state=True) + return m.fs.state[0] + + +def initialize(m): + assert degrees_of_freedom(m) == 0 + m.fs.state.initialize() + assert degrees_of_freedom(m) == 0 + + def solve(m): solver = SolverFactory('ipopt') solver.solve(m, tee=True) + +def test_constrain(): + m = flowsheet() + sb = build_state(m) + # fix flow_mol directly + c = sb.constrain("flow_mol", 1) + assert c == sb.flow_mol + assert c.value == 1 + assert c.is_fixed() + + # add a constraint for flow_mass + c = sb.constrain("flow_mass", 1) + assert type(c) == ScalarConstraint + assert c in sb.component_data_objects(Constraint) + assert getattr(sb.constraints, "flow_mass") == c + + def test_state_definition(): m = flowsheet() - m.fs.sb = m.fs.properties.build_state_block(m.fs.time) - sb = m.fs.sb[0] + sb = build_state(m) sb.constrain("enth_mass", 1878.87) sb.constrain("pressure", 101325) sb.constrain("flow_mass", 1) - m.fs.sb.initialize() + initialize(m) solve(m) assert value(sb.temperature) == approx(273.5809) def test_state_definition_temp(): m = flowsheet() - m.fs.sb = m.fs.properties.build_state_block(m.fs.time) - sb = m.fs.sb[0] - sb.constrain("smooth_temperature", 273.5809) + sb = build_state(m) + sb.constrain("temperature", 273.5809) sb.constrain("pressure", 101325) sb.constrain("flow_mass", 1) - m.fs.sb.initialize() + initialize(m) solve(m) assert value(sb.enth_mass) == approx(1878.712) @@ -48,12 +77,13 @@ def test_initialise_compressor(): m.fs.compressor = Compressor(property_package=m.fs.properties) inlet = m.fs.compressor.control_volume.properties_in[0] outlet = m.fs.compressor.control_volume.properties_out[0] - inlet.constrain("smooth_temperature", 273.5809) + inlet.constrain("temperature", 273.5809) inlet.constrain("pressure", 101325) inlet.constrain("flow_mass", 1) + m.fs.compressor.efficiency_isentropic.fix(1) m.fs.compressor.deltaP.fix(100000) - m.fs.compressor.initialize() + m.fs.compressor.initialize(outlvl=1) solve(m) - assert value(outlet.temperature) == approx(393.5689573) + assert value(outlet.temperature) == approx(273.58047830928655) diff --git a/property_packages/modular/builder/common_parsers.py b/property_packages/modular/builder/common_parsers.py index 2c99e3b..6360323 100644 --- a/property_packages/modular/builder/common_parsers.py +++ b/property_packages/modular/builder/common_parsers.py @@ -3,7 +3,7 @@ from .base_parser import BuildBase from compounds.Compound import Compound from pyomo.environ import units as pyunits -from idaes.models.properties.modular_properties.state_definitions import FTPx +from idaes.models.properties.modular_properties.state_definitions import FTPx, FPhx from idaes.models.properties.modular_properties.phase_equil.bubble_dew import (LogBubbleDew) from idaes.core import LiquidPhase, VaporPhase, Component, PhaseType as PT from idaes.models.properties.modular_properties.phase_equil import (SmoothVLE) @@ -246,30 +246,12 @@ def serialise(compounds: List[Compound], valid_states: List[States]) -> Dict[str min_critical_temperature = min([compound["CriticalTemperature"].value for compound in compounds]) # TODO: Refactor this logic, need a more versatile approach - if (min_critical_temperature > 500): - return { - "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), - "temperature": (min_melting_point, 300, 500, pyunits.K), - "pressure": (5e4, 1e5, 1e6, pyunits.Pa), - } - elif (min_critical_temperature > 300): - return { - "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), - "temperature": (min_melting_point, 200, 400, pyunits.K), - "pressure": (5e4, 1e5, 1e6, pyunits.Pa), - } - elif (min_critical_temperature > 120): - return { - "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), - "temperature": (min_melting_point, 150, 500, pyunits.K), - "pressure": (5e4, 1e5, 1e6, pyunits.Pa), - } - else: - return { - "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), - "temperature": (min_melting_point, 150, 350, pyunits.K), - "pressure": (5e4, 1e5, 1e6, pyunits.Pa), - } + return { + "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), + "temperature": (max(min_melting_point-50,1), 300, 3000, pyunits.K), + "pressure": (5e4, 1e5, 1e6, pyunits.Pa), + } + class state_definition_parser(BuildBase): @staticmethod diff --git a/property_packages/modular/modular_extended.py b/property_packages/modular/modular_extended.py index cff816b..298a146 100644 --- a/property_packages/modular/modular_extended.py +++ b/property_packages/modular/modular_extended.py @@ -1,115 +1,552 @@ -# extends the IDAES helmholtz property package to include additional properties and methods. -from idaes.models.properties.general_helmholtz.helmholtz_state import HelmholtzStateBlockData, _StateBlock -from idaes.models.properties.general_helmholtz.helmholtz_functions import HelmholtzParameterBlockData -from idaes.models.properties.modular_properties.base.generic_property import GenericParameterBlock, _GenericStateBlock, GenericParameterData, GenericStateBlockData from idaes.core import declare_process_block_class +from idaes.core.util.exceptions import ConfigurationError +from idaes.models.properties.modular_properties.base.generic_property import ( + _GenericStateBlock, + GenericParameterData, + GenericStateBlockData, +) +import idaes.models.properties.modular_properties.base.utility as utility + +from property_packages.base.state_block_constraints import StateBlockConstraints from property_packages.utils.add_extra_expressions import add_extra_expressions -from pyomo.environ import Constraint, Block -from pyomo.core.base.expression import Expression, ScalarExpression, _GeneralExpressionData, ExpressionData -from pyomo.core.base.var import IndexedVar, ScalarVar, Var, _GeneralVarData,VarData +from property_packages.utils.fix_state_vars import fix_state_vars + + +# TODO: remove these imports once https://github.com/IDAES/idaes-pse/pull/1554 is resolved +from pyomo.environ import ( + Block, + check_optimal_termination, + Constraint, + exp, + Expression, + log, + Set, + Param, + value, + Var, + units as pyunits, + Reference, +) +from pyomo.common.config import ConfigBlock, ConfigDict, ConfigValue, In, Bool +from pyomo.util.calc_var_value import calculate_variable_from_constraint + +# Import IDAES cores +from idaes.core import ( + declare_process_block_class, + PhysicalParameterBlock, + StateBlockData, + StateBlock, + MaterialFlowBasis, + ElectrolytePropertySet, +) +from idaes.core.base.components import Component, __all_components__ +from idaes.core.base.phases import ( + Phase, + AqueousPhase, + LiquidPhase, + VaporPhase, + __all_phases__, +) +from idaes.core.util.initialization import ( + revert_state_vars, + solve_indexed_blocks, +) +from idaes.core.util.model_statistics import ( + degrees_of_freedom, + number_activated_constraints, +) +from idaes.core.util.exceptions import ( + BurntToast, + ConfigurationError, + PropertyPackageError, + PropertyNotSupportedError, + InitializationError, +) +from idaes.core.util.misc import add_object_reference +from idaes.core.solvers import get_solver import idaes.logger as idaeslog +import idaes.core.util.scaling as iscale +from idaes.core.initialization.initializer_base import InitializerBase -# NOTE: -# THis only works for FTPx formulation right now. +from idaes.models.properties.modular_properties.base.generic_reaction import ( + equil_rxn_config, +) +from idaes.models.properties.modular_properties.base.utility import ( + get_method, + get_phase_method, + GenericPropertyPackageError, + StateIndex, + identify_VL_component_list, + estimate_Tbub, + estimate_Tdew, + estimate_Pbub, + estimate_Pdew, +) +from idaes.models.properties.modular_properties.phase_equil.bubble_dew import ( + LogBubbleDew, +) +from idaes.models.properties.modular_properties.phase_equil.henry import HenryType +from idaes.models.properties.modular_properties.base.generic_property import ( + _initialize_critical_props, + _init_Tbub, + _init_Tdew, + _init_Pbub, + _init_Pdew, +) + + +# increase max iterations for estimating values for bubble, dew, and +# critical point initialization (ie. temperature_bubble, temperature_dew) +utility.MAX_ITER = 1000 class _ExtendedGenericStateBlock(_GenericStateBlock): - """ - This class contains methods which should be applied to Property Blocks as a - whole, rather than individual elements of indexed Property Blocks. - """ - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) # Missing argument - - def initialize(self, *args, **kwargs): - print("GETTING STATE") - hold_state = kwargs.pop("hold_state", False) - for i, v in self.items(): - print(f"State block {i}") - print(f"block data: {v}") - v.constraints.deactivate() - print("ACTIVATING CONSTRAINTS") - res = super().initialize(*args, **kwargs) - flags = {} - print("STATE BLOCKS") - print(self.items()) - for i, v in self.items(): - print(f"STATE {i}") - print(f"BLOCK DATA {v}") - v.constraints.activate() - flags[i] = {} - if hold_state: - # Fix the required state variables for zero degrees of freedom, and return a dictionary of the flags. - if not hasattr(v.constraints,"flow_mass") and not v.flow_mol.is_fixed(): - # We need to fix the flow_mol variable - flags[i]["flow_mol"] = True - v.flow_mol.fix() - avaliable_constraints = ["enth_mass","temperature","entr_mass","entr_mol","smooth_temperature","vapor_frac"] - if not v.enth_mol.is_fixed(): - # check if any of the constraints exist - found_constraint = False - for constraint in avaliable_constraints: - if hasattr(v.constraints,constraint): - # we don't need to fix the variable - # but we need to remove this from the list of constraints (it can't be used to define pressure) - avaliable_constraints.remove(constraint) - found_constraint = True - break - if not found_constraint: - # we need to fix the variable - flags[i]["enth_mol"] = True - v.enth_mol.fix() - if not v.pressure.is_fixed(): - # check if any of the constraints exist - found_constraint = False - for constraint in avaliable_constraints: - if hasattr(v.constraints,constraint): - # we don't need to fix the variable - # but we need to remove this from the list of constraints (it can't be used to define pressure) - avaliable_constraints.remove(constraint) - found_constraint = True - break - if not found_constraint: - # we need to fix the variable - flags[i]["pressure"] = True - v.pressure.fix() - print("WE DONE") - return flags - - def release_state(self, flags, outlvl=idaeslog.NOTSET): - print("RELEASING STATE") - for i, v in self.items(): - for key in flags[i]: - getattr(v,key).unfix() - -@declare_process_block_class("GenericExtendedStateBlock", block_class=_ExtendedGenericStateBlock) -class GenericExtendedStateBlockData(GenericStateBlockData): + def initialize(blk, *args, **kwargs): + flag_dict = fix_state_vars(blk, kwargs.get("state_args", None)) - def build(self, *args): - super().build(*args) - # Add expressions for smooth_temperature, enthalpy in terms of mass, etc. - add_extra_expressions(self) - # Add a block for constraints, so we can disable or enable them in bulk - self.constraints = Block() + # Set state_vars_fixed to True to avoid fixing state variables + # during the initialize method, since this would overdefine + # the block if we are using constraints + kwargs["state_vars_fixed"] = True + + # TODO: replace this with + # super().initialize(*args, **kwargs) + # once https://github.com/IDAES/idaes-pse/pull/1554 has + # been resolved and released (allows using constraints to + # define state variables during initialization) + blk._custom_super_initialize(*args, **kwargs) + + if kwargs.get("hold_state") is True: + return flag_dict + else: + blk.release_state(flag_dict) - def constrain(self,name:str,value:float): - # Value must be a float. TODO: Handle unit conversion. - var = getattr(self,name) - if type(var) == ScalarExpression: - self.constraints.add_component(name, Constraint(expr=var == value)) - elif type(var) in (ScalarVar, _GeneralVarData, VarData): - var.fix(value) - elif type(var) in ( _GeneralExpressionData, ExpressionData) : - # allowed, but we don't need to fix it (eg. mole_frac_comp in helmholtz) - print(f"Variable {self} {name} is an Expression: {type(var)}") - pass + + def _custom_super_initialize( + blk, + state_args=None, + state_vars_fixed=False, + hold_state=False, + outlvl=idaeslog.NOTSET, + solver=None, + optarg=None, + ): + """ + Initialization routine for property package. + Keyword Arguments: + state_args : a dict of initial values for the state variables + defined by the property package. + outlvl : sets output level of initialization routine + optarg : solver options dictionary object (default=None, use + default solver options) + state_vars_fixed: Flag to denote if state vars have already been + fixed. + - True - states have already been fixed by the + control volume 1D. Control volume 0D + does not fix the state vars, so will + be False if this state block is used + with 0D blocks. + - False - states have not been fixed. The state + block will deal with fixing/unfixing. + solver : str indicating which solver to use during + initialization (default = None, use default solver) + hold_state : flag indicating whether the initialization routine + should unfix any state variables fixed during + initialization (default=False). + - True - states variables are not unfixed, and + a dict of returned containing flags for + which states were fixed during + initialization. + - False - state variables are unfixed after + initialization by calling the + release_state method + Returns: + If hold_states is True, returns a dict containing flags for + which states were fixed during initialization. + """ + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties") + + init_log.info("Starting initialization") + + res = None + + for k in blk.values(): + # Deactivate the constraints specific for outlet block i.e. + # when defined state is False + if k.config.defined_state is False: + try: + k.sum_mole_frac_out.deactivate() + except AttributeError: + pass + + if hasattr(k, "inherent_equilibrium_constraint") and ( + not k.params._electrolyte + or k.params.config.state_components == StateIndex.true + ): + k.inherent_equilibrium_constraint.deactivate() + + # Fix state variables if not already fixed + if state_vars_fixed is False: + flag_dict = fix_state_vars(blk, state_args) + # Confirm DoF for sanity + for k in blk.values(): + if k.always_flash: + # If not always flash, DoF is probably less than zero + # We will handle this elsewhere + dof = degrees_of_freedom(k) + if dof != 0: + raise BurntToast( + "Degrees of freedom were not zero [{}] " + "after trying to fix state variables. " + "Something broke in the generic property " + "package code - please inform the IDAES " + "developers.".format(dof) + ) else: - raise Exception(f"Variable {self} {name} is not a Var or Expression: {type(var)}") + # When state vars are fixed, check that DoF is 0 + for k in blk.values(): + if degrees_of_freedom(k) != 0: + # PYLINT-TODO + # pylint: disable-next=broad-exception-raised + raise Exception( + "State vars fixed but degrees of " + "freedom for state block is not zero " + "during initialization." + ) + + # Create solver + opt = get_solver(solver, optarg) + + # --------------------------------------------------------------------- + # If present, initialize bubble, dew , and critical point calculations + for k in blk.values(): + T_units = k.params.get_metadata().default_units.TEMPERATURE + + # List of bubble and dew point constraints + cons_list = [ + "eq_pressure_dew", + "eq_pressure_bubble", + "eq_temperature_dew", + "eq_temperature_bubble", + "eq_mole_frac_tbub", + "eq_mole_frac_tdew", + "eq_mole_frac_pbub", + "eq_mole_frac_pdew", + "log_mole_frac_tbub_eqn", + "log_mole_frac_tdew_eqn", + "log_mole_frac_pbub_eqn", + "log_mole_frac_pdew_eqn", + "mole_frac_comp_eq", + "log_mole_frac_comp_eqn", + ] + + # Critical point + with k.lock_attribute_creation_context(): + # Only need to look for one, as it is all-or-nothing + if hasattr(k, "pressure_crit"): + # Initialize critical point properties + _initialize_critical_props(k) + # Add critical point constraints to cons_list + cons_list += k.list_critical_property_constraint_names() + + # Bubble temperature initialization + if hasattr(k, "_mole_frac_tbub"): + _init_Tbub(k, T_units) + + # Dew temperature initialization + if hasattr(k, "_mole_frac_tdew"): + _init_Tdew(k, T_units) + + # Bubble pressure initialization + if hasattr(k, "_mole_frac_pbub"): + _init_Pbub(k) + + # Dew pressure initialization + if hasattr(k, "_mole_frac_pdew"): + _init_Pdew(k) + + # Solve bubble, dew, and critical point constraints + for c in k.component_objects(Constraint): + # Deactivate all constraints not associated with bubble, dew, + # or critical points + if c.local_name not in cons_list: + c.deactivate() + + # If StateBlock has active constraints (i.e. has bubble, dew, or critical + # point calculations), solve the block to converge these + n_cons = 0 + dof = 0 + for k in blk.values(): + n_cons += number_activated_constraints(k) + dof += degrees_of_freedom(k) + if n_cons > 0: + if dof > 0: + raise InitializationError( + f"{blk.name} Unexpected degrees of freedom during " + f"initialization at bubble, dew, and critical point step: {dof}." + ) + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solve_indexed_blocks(opt, [blk], tee=slc.tee) + init_log.info( + "Bubble, dew, and critical point initialization: {}.".format( + idaeslog.condition(res) + ) + ) + # --------------------------------------------------------------------- + # Calculate _teq if required + # Using iterator k outside of for loop - this should be OK as we just need + # a valid StateBlockData an assume they are all the same. + if k.params.config.phases_in_equilibrium is not None and ( + not k.config.defined_state or k.always_flash + ): + for k in blk.values(): + for pp in k.params._pe_pairs: + k.params.config.phase_equilibrium_state[pp].calculate_teq(k, pp) + + init_log.info("Equilibrium temperature initialization completed.") + + # --------------------------------------------------------------------- + # Initialize flow rates and compositions + for k in blk.values(): + + k.params.config.state_definition.state_initialization(k) + + if k.params._electrolyte: + if k.params.config.state_components == StateIndex.true: + # First calculate initial values for apparent species flows + for p, j in k.params.apparent_phase_component_set: + calculate_variable_from_constraint( + k.flow_mol_phase_comp_apparent[p, j], + k.true_to_appr_species[p, j], + ) + # Need to calculate all flows before doing mole fractions + for p, j in k.params.apparent_phase_component_set: + sum_flow = sum( + k.flow_mol_phase_comp_apparent[p, jj] + for jj in k.params.apparent_species_set + if (p, jj) in k.params.apparent_phase_component_set + ) + if value(sum_flow) == 0: + x = 1 + else: + x = value(k.flow_mol_phase_comp_apparent[p, j] / sum_flow) + lb = k.mole_frac_phase_comp_apparent[p, j].lb + if lb is not None and x <= lb: + k.mole_frac_phase_comp_apparent[p, j].set_value(lb) + else: + k.mole_frac_phase_comp_apparent[p, j].set_value(x) + elif k.params.config.state_components == StateIndex.apparent: + # First calculate initial values for true species flows + for p, j in k.params.true_phase_component_set: + calculate_variable_from_constraint( + k.flow_mol_phase_comp_true[p, j], + k.appr_to_true_species[p, j], + ) + # Need to calculate all flows before doing mole fractions + for p, j in k.params.true_phase_component_set: + sum_flow = sum( + k.flow_mol_phase_comp_true[p, jj] + for jj in k.params.true_species_set + if (p, jj) in k.params.true_phase_component_set + ) + if value(sum_flow) == 0: + x = 1 + else: + x = value(k.flow_mol_phase_comp_true[p, j] / sum_flow) + lb = k.mole_frac_phase_comp_true[p, j].lb + if lb is not None and x <= lb: + k.mole_frac_phase_comp_true[p, j].set_value(lb) + else: + k.mole_frac_phase_comp_true[p, j].set_value(x) + + # If state block has phase equilibrium, use the average of all + # _teq's as an initial guess for T + if ( + k.params.config.phases_in_equilibrium is not None + and isinstance(k.temperature, Var) + and not k.temperature.fixed + ): + k.temperature.value = value( + sum(k._teq[i] for i in k.params._pe_pairs) / len(k.params._pe_pairs) + ) + + if outlvl > 0: # TODO: Update to use logger Enum + init_log.info("State variable initialization completed.") + + # --------------------------------------------------------------------- + n_cons = 0 + dof = 0 + skip = False + Tfix = {} # In enth based state defs, need to also fix T until later + for k, b in blk.items(): + if b.params.config.phase_equilibrium_state is not None and ( + not b.config.defined_state or b.always_flash + ): + if not b.temperature.fixed: + b.temperature.fix() + Tfix[k] = True + for c in b.component_objects(Constraint): + # Activate common constraints + if c.local_name in ( + "total_flow_balance", + "component_flow_balances", + "sum_mole_frac", + "phase_fraction_constraint", + "mole_frac_phase_comp_eq", + "mole_frac_comp_eq", + ): + c.activate() + if c.local_name == "log_mole_frac_phase_comp_eqn": + c.activate() + for p, j in b.params._phase_component_set: + calculate_variable_from_constraint( + b.log_mole_frac_phase_comp[p, j], + b.log_mole_frac_phase_comp_eqn[p, j], + ) + elif c.local_name == "equilibrium_constraint": + # 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(): + c.activate() + elif getattr(c, "defining_state_var", False): + # Allow using a constraint to define a state var + # (rather than fixing it directly) by setting + # c.defining_state_var = True. This can be used in + # conjunction with setting state_vars_fixed = True + c.activate() + + for pp in b.params._pe_pairs: + # Activate formulation specific constraints + b.params.config.phase_equilibrium_state[ + pp + ].phase_equil_initialization(b, pp) + + n_cons += number_activated_constraints(b) + dof += degrees_of_freedom(b) + if degrees_of_freedom(b) < 0: + # Skip solve if DoF < 0 - this is probably due to a + # phase-component flow state with flash + skip = True + + if n_cons > 0 and not skip: + if dof > 0: + raise InitializationError( + f"{blk.name} Unexpected degrees of freedom during " + f"initialization at phase equilibrium step: {dof}." + ) + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solve_indexed_blocks(opt, [blk], tee=slc.tee) + init_log.info( + "Phase equilibrium initialization: {}.".format(idaeslog.condition(res)) + ) + + # --------------------------------------------------------------------- + # Initialize other properties + for k, b in blk.items(): + for c in b.component_objects(Constraint): + # Activate all constraints except flagged do_not_initialize + if c.local_name not in ( + b.params.config.state_definition.do_not_initialize + ): + c.activate() + if k in Tfix: + b.temperature.unfix() + + # Initialize log-form variables + log_form_vars = [ + "act_phase_comp", + "act_phase_comp_apparent", + "act_phase_comp_true", + "conc_mol_phase_comp", + "conc_mol_phase_comp_apparent", + "conc_mol_phase_comp_true", + "mass_frac_phase_comp", + "mass_frac_phase_comp_apparent", + "mass_frac_phase_comp_true", + "molality_phase_comp", + "molality_phase_comp_apparent", + "molality_phase_comp_true", + "mole_frac_comp", # Might have already been initialized + "mole_frac_phase_comp", # Might have already been initialized + "mole_frac_phase_comp_apparent", + "mole_frac_phase_comp_true", + "pressure_phase_comp", + "pressure_phase_comp_apparent", + "pressure_phase_comp_true", + ] + + for prop in log_form_vars: + if b.is_property_constructed("log_" + prop): + comp = getattr(b, prop) + lcomp = getattr(b, "log_" + prop) + for k2, v in lcomp.items(): + c = value(comp[k2]) + if c <= 0: + c = 1e-8 + lc = log(c) + v.set_value(value(lc)) + + n_cons = 0 + dof = 0 + skip = False + for k in blk.values(): + if degrees_of_freedom(k) < 0: + # Skip solve if DoF < 0 - this is probably due to a + # phase-component flow state with flash + skip = True + n_cons += number_activated_constraints(k) + dof += degrees_of_freedom(k) + if n_cons > 0 and not skip: + if dof > 0: + raise InitializationError( + f"{blk.name} Unexpected degrees of freedom during " + f"initialization at property initialization step: {dof}." + ) + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solve_indexed_blocks(opt, [blk], tee=slc.tee) + init_log.info( + "Property initialization: {}.".format(idaeslog.condition(res)) + ) + + # --------------------------------------------------------------------- + # Return constraints to initial state + for k in blk.values(): + for c in k.component_objects(Constraint): + if c.local_name in (k.params.config.state_definition.do_not_initialize): + c.activate() + + if res is not None and not check_optimal_termination(res): + raise InitializationError( + f"{blk.name} failed to initialize successfully. Please check " + f"the output logs for more information." + ) + + if state_vars_fixed is False: + if hold_state is True: + return flag_dict + else: + blk.release_state(flag_dict) + + init_log.info( + "Property package initialization: {}.".format(idaeslog.condition(res)) + ) + + +@declare_process_block_class( + "GenericExtendedStateBlock", block_class=_ExtendedGenericStateBlock +) +class GenericExtendedStateBlockData(GenericStateBlockData, StateBlockConstraints): + + def build(self, *args): + GenericStateBlockData.build(self, *args) + StateBlockConstraints.build(self, *args) + @declare_process_block_class("GenericExtendedParameterBlock") class GenericExtendedParameterData(GenericParameterData): def build(self): super().build() - # set that we should use the extended state block - self._state_block_class = GenericExtendedStateBlock # type: ignore because it'll get created with declare_process_block_class + self._state_block_class = GenericExtendedStateBlock # noqa: F821 diff --git a/property_packages/modular/template_builder.py b/property_packages/modular/template_builder.py index 7e1c091..5aec4dc 100644 --- a/property_packages/modular/template_builder.py +++ b/property_packages/modular/template_builder.py @@ -1,5 +1,6 @@ from idaes.models.properties.modular_properties.base.generic_property import GenericParameterBlock +from .modular_extended import GenericExtendedParameterBlock from compounds.CompoundDB import get_compound from .templates.templates import PropertyPackage from property_packages.types import States @@ -44,4 +45,4 @@ def build_config(property_package_name, compound_names: List[str], valid_states: new_template[key] = obj.serialise(compounds, valid_states) # Building property package and returning - return GenericParameterBlock(**new_template) + return GenericExtendedParameterBlock(**new_template) diff --git a/property_packages/tests/test_asu_pr.py b/property_packages/modular/tests/test_asu_pr.py similarity index 96% rename from property_packages/tests/test_asu_pr.py rename to property_packages/modular/tests/test_asu_pr.py index 5088c2b..c46e0c3 100644 --- a/property_packages/tests/test_asu_pr.py +++ b/property_packages/modular/tests/test_asu_pr.py @@ -23,7 +23,7 @@ from idaes.models.properties.modular_properties.examples.ASU_PR import configuration from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available from idaes.models.properties.modular_properties.base.generic_property import GenericParameterBlock -from ..build_package import build_package +from property_packages.build_package import build_package solver = get_solver("ipopt") @@ -76,7 +76,7 @@ def test_build(self): model.params.config.state_bounds, { "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), - "temperature": (54.361, 150, 500, pyunits.K), + "temperature": (4.361, 300, 3000, pyunits.K), "pressure": (5e4, 1e5, 1e6, pyunits.Pa), }, item_callback=_as_quantity, @@ -138,9 +138,9 @@ def test_build(self, model): assert model.props[1].pressure.lb == 5e4 assert isinstance(model.props[1].temperature, Var) - assert value(model.props[1].temperature) == 150 - assert model.props[1].temperature.ub == 500 - assert model.props[1].temperature.lb == 54.361 + assert value(model.props[1].temperature) == 300 + assert model.props[1].temperature.ub == 3000 + assert_approx(model.props[1].temperature.lb, 4.361, 0.1) assert isinstance(model.props[1].mole_frac_comp, Var) assert len(model.props[1].mole_frac_comp) == 3 diff --git a/property_packages/modular/tests/test_bt_pr.py b/property_packages/modular/tests/test_bt_pr.py new file mode 100644 index 0000000..5ccfaf6 --- /dev/null +++ b/property_packages/modular/tests/test_bt_pr.py @@ -0,0 +1,457 @@ +# Build and solve a state block. +from property_packages.build_package import build_package +from pytest import approx + +# Import objects from pyomo package +from pyomo.environ import ConcreteModel, value + +# Import the main FlowsheetBlock from IDAES. The flowsheet block will contain the unit model +from idaes.core import FlowsheetBlock +from idaes.core.solvers import get_solver + +from pyomo.util.check_units import assert_units_consistent + +from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available +from pyomo.environ import check_optimal_termination, ConcreteModel, Objective +import idaes.core.util.scaling as iscale +from idaes.core.util.model_statistics import degrees_of_freedom +solver = get_solver(solver="ipopt") + +import idaes.logger as idaeslog +SOUT = idaeslog.INFO + +""" +Test Suite Sourced From IDAES + +URL: https://github.com/IDAES/idaes-pse/blob/41bb3c9728ea227fa8bb0fa1bf35b8deec467783/idaes/models/ +properties/modular_properties/examples/tests/test_BT_PR.py +""" + +""" +--------------------------------------------- +Goal for absolute error margins +--------------------------------------------- +- temp, pressure, specific volume : 0.5% +- compositional things : 2% (mole frac, flow, etc) +- thermodynamic values: 5% (entr, enth) +--------------------------------------------- +""" + +""" +Helper function to calculate the absolute error margin of a +value and asserts whether or not value lies within range +- Accepts: percent_error (as whole percent) +- Returns: None +""" +def assert_approx(value, expected_value, error_margin): + percent_error = error_margin / 100 + tolerance = abs(percent_error * expected_value) + assert approx(value, abs=tolerance) == expected_value + +def get_m(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False,time_set=[1], ) + m.fs.props = build_package("peng-robinson", ["benzene", "toluene"], ["Liq", "Vap"]) + m.fs.state = m.fs.props.build_state_block([1], defined_state=True) + iscale.calculate_scaling_factors(m.fs.props) + iscale.calculate_scaling_factors(m.fs.state[1]) + return m + +def test_T_sweep(): + m = get_m() + sb = m.fs.state[1] + + assert_units_consistent(m) + + m.fs.obj = Objective(expr=(sb.temperature - 510) ** 2) + sb.temperature.setub(600) + + # Tests a variety of pressures, and makes sure that the benzene-toluene mixture + # can move from liquid to vapor phase at each pressure + for logP in [9.5, 10, 10.5, 11, 11.5, 12]: + m.fs.obj.deactivate() + + sb.flow_mol.fix(100) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + sb.pressure.fix(10 ** (0.5 * logP)) + + assert degrees_of_freedom(sb) == 1 + + m.fs.state.initialize() + sb.temperature.fix(300) + + assert degrees_of_freedom(sb) == 0 + + results = solver.solve(m) + assert sb.flow_mol_phase["Vap"].value <= 1e-2 + + sb.temperature.unfix() + assert degrees_of_freedom(sb) == 1 + m.fs.obj.activate() + + results = solver.solve(m) + + assert check_optimal_termination(results) + assert sb.flow_mol_phase["Liq"].value <= 1e-2 + + +def test_P_sweep(): + m = get_m() + sb = m.fs.state[1] + + # Tests initiialization at a variety of temperatures. + for T in range(370, 500, 25): + sb.flow_mol.fix(100) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + sb.pressure.fix(1e5) + sb.temperature.fix(T) + m.fs.state.initialize() + + results = solver.solve(m) + + assert check_optimal_termination(results) + + while sb.pressure.value <= 1e6: + + results = solver.solve(m) + assert check_optimal_termination(results) + + sb.pressure.value = sb.pressure.value + 1e5 + +def test_T350_P1_x5(): + m = get_m() + sb = m.fs.state[1] + + sb.flow_mol.fix(100) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + sb.temperature.fix(350) + sb.pressure.fix(1e5) + + # Trigger build of enthalpy and entropy + sb.enth_mol_phase + sb.entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert_approx(value(sb._teq[("Vap", "Liq")]), 365, 0.5) + assert_approx(value(sb.compress_fact_phase["Liq"]), 0.0035346, 0.5) + assert_approx(value(sb.compress_fact_phase["Vap"]), 0.966749, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "benzene"]), 0.894676, 2) # Investigate + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "toluene"]), 0.347566, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "benzene"]), 0.971072, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "toluene"]), 0.959791, 0.5) + + assert_approx(value(sb.mole_frac_phase_comp["Liq", "benzene"]), 0.5, 2) + assert_approx(value(sb.mole_frac_phase_comp["Liq", "toluene"]), 0.5, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "benzene"]), 0.70584, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "toluene"]), 0.29416, 2) + + assert_approx(value(sb.enth_mol_phase["Liq"]), 38942.8, 5) + assert_approx(value(sb.enth_mol_phase["Vap"]), 78048.7, 5) + assert_approx(value(sb.entr_mol_phase["Liq"]), -361.794, 5) + assert_approx(value(sb.entr_mol_phase["Vap"]), -264.0181, 5) + +def test_T350_P5_x5(): + m = get_m() + sb = m.fs.state[1] + + sb.flow_mol.fix(100) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + sb.temperature.fix(350) + sb.pressure.fix(5e5) + + # Trigger build of enthalpy and entropy + sb.enth_mol_phase + sb.entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert_approx(value(sb._teq[("Vap", "Liq")]), 431.47, 0.5) + assert_approx(value(sb.compress_fact_phase["Liq"]), 0.01766, 0.5) + assert_approx(value(sb.compress_fact_phase["Vap"]), 0.80245, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "benzene"]), 0.181229, 1.5) # Investigate + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "toluene"]), 0.070601, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "benzene"]), 0.856523, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "toluene"]), 0.799237, 0.5) + + assert_approx(value(sb.mole_frac_phase_comp["Liq", "benzene"]), 0.5, 2) + assert_approx(value(sb.mole_frac_phase_comp["Liq", "toluene"]), 0.5, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "benzene"]), 0.65415, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "toluene"]), 0.34585, 2) + + assert_approx(value(sb.enth_mol_phase["Liq"]), 38966.9, 5) + assert_approx(value(sb.enth_mol_phase["Vap"]), 75150.7, 5) + assert_approx(value(sb.entr_mol_phase["Liq"]), -361.8433, 5) + assert_approx(value(sb.entr_mol_phase["Vap"]), -281.9703, 5) + +def test_T450_P1_x5(): + m = get_m() + sb = m.fs.state[1] + + sb.flow_mol.fix(100) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + sb.temperature.fix(450) + sb.pressure.fix(1e5) + + # Trigger build of enthalpy and entropy + sb.enth_mol_phase + sb.entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert_approx(value(sb._teq[("Vap", "Liq")]), 371.4, 0.5) + assert_approx(value(sb.compress_fact_phase["Liq"]), 0.0033583, 0.5) + assert_approx(value(sb.compress_fact_phase["Vap"]), 0.9821368, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "benzene"]), 8.069323, 1) # Investigate + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "toluene"]), 4.304955, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "benzene"]), 0.985365, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "toluene"]), 0.979457, 0.5) + + assert_approx(value(sb.mole_frac_phase_comp["Liq", "benzene"]), 0.29861, 2) + assert_approx(value(sb.mole_frac_phase_comp["Liq", "toluene"]), 0.70139, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "benzene"]), 0.5, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "toluene"]), 0.5, 2) + + assert_approx(value(sb.enth_mol_phase["Liq"]), 49441.2, 5) + assert_approx(value(sb.enth_mol_phase["Vap"]), 84175.1, 5) + assert_approx(value(sb.entr_mol_phase["Liq"]), -328.766, 5) + assert_approx(value(sb.entr_mol_phase["Vap"]), -241.622, 5) + +def test_T450_P5_x5(): + m = get_m() + sb = m.fs.state[1] + + sb.flow_mol.fix(100) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + sb.temperature.fix(450) + sb.pressure.fix(5e5) + + # Trigger build of enthalpy and entropy + sb.enth_mol_phase + sb.entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert_approx(value(sb._teq[("Vap", "Liq")]), 436.93, 0.5) + assert_approx(value(sb.compress_fact_phase["Liq"]), 0.0166181, 0.5) + assert_approx(value(sb.compress_fact_phase["Vap"]), 0.9053766, 0.5) + + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "benzene"]), 1.63308, 1) # investigate + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "toluene"]), 0.873213, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "benzene"]), 0.927534, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "toluene"]), 0.898324, 0.5) + + assert_approx(value(sb.mole_frac_phase_comp["Liq", "benzene"]), 0.3488737, 2) + assert_approx(value(sb.mole_frac_phase_comp["Liq", "toluene"]), 0.6511263, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "benzene"]), 0.5, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "toluene"]), 0.5, 2) + + assert_approx(value(sb.enth_mol_phase["Liq"]), 51095.2, 5) + assert_approx(value(sb.enth_mol_phase["Vap"]), 83362.3, 5) + assert_approx(value(sb.entr_mol_phase["Liq"]), -326.299, 5) + assert_approx(value(sb.entr_mol_phase["Vap"]), -256.198, 5) + +def test_T368_P1_x5(): + m = get_m() + sb = m.fs.state[1] + + sb.flow_mol.fix(100) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + sb.temperature.fix(368) + sb.pressure.fix(1e5) + + # Trigger build of enthalpy and entropy + sb.enth_mol_phase + sb.entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert_approx(value(sb._teq[("Vap", "Liq")]), 368, 0.5) + assert_approx(value(sb.compress_fact_phase["Liq"]), 0.003504, 0.5) + assert_approx(value(sb.compress_fact_phase["Vap"]), 0.97, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "benzene"]), 1.492049, 1.5) # Investigate + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "toluene"]), 0.621563, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "benzene"]), 0.97469, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "toluene"]), 0.964642, 0.5) + + assert_approx(value(sb.mole_frac_phase_comp["Liq", "benzene"]), 0.4012128, 2) + assert_approx(value(sb.mole_frac_phase_comp["Liq", "toluene"]), 0.5987872, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "benzene"]), 0.6141738, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "toluene"]), 0.3858262, 2) + + assert_approx(value(sb.enth_mol_phase["Liq"]), 38235.1, 5) + assert_approx(value(sb.enth_mol_phase["Vap"]), 77155.4, 5) + assert_approx(value(sb.entr_mol_phase["Liq"]), -359.256, 5) + assert_approx(value(sb.entr_mol_phase["Vap"]), -262.348, 5) + +def test_T376_P1_x2(): + m = get_m() + sb = m.fs.state[1] + + sb.flow_mol.fix(100) + sb.mole_frac_comp["benzene"].fix(0.2) + sb.mole_frac_comp["toluene"].fix(0.8) + sb.temperature.fix(376) + sb.pressure.fix(1e5) + + # Trigger build of enthalpy and entropy + sb.enth_mol_phase + sb.entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert_approx(value(sb._teq[("Vap", "Liq")]), 376, 0.5) + assert_approx(value(sb.compress_fact_phase["Liq"]), 0.00361333, 0.5) + assert_approx(value(sb.compress_fact_phase["Vap"]), 0.968749, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "benzene"]), 1.8394188, 1.5) # Investigate + assert_approx(value(sb.fug_coeff_phase_comp["Liq", "toluene"]), 0.7871415, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "benzene"]), 0.9763608, 0.5) + assert_approx(value(sb.fug_coeff_phase_comp["Vap", "toluene"]), 0.9663611, 0.5) + + assert_approx(value(sb.mole_frac_phase_comp["Liq", "benzene"]), 0.17342, 2) + assert_approx(value(sb.mole_frac_phase_comp["Liq", "toluene"]), 0.82658, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "benzene"]), 0.3267155, 2) + assert_approx(value(sb.mole_frac_phase_comp["Vap", "toluene"]), 0.6732845, 2) + + assert_approx(value(sb.enth_mol_phase["Liq"]), 31535.8, 5) + assert_approx(value(sb.enth_mol_phase["Vap"]), 69175.3, 5) + assert_approx(value(sb.entr_mol_phase["Liq"]), -369.033, 5) + assert_approx(value(sb.entr_mol_phase["Vap"]), -273.513, 5) + +def test_basic_scaling(): + m = get_m() + sb = m.fs.state[1] + assert len(sb.scaling_factor) == 26 # IDAES version has 23, but I think we have some more because of enth_mol + assert sb.scaling_factor[sb.flow_mol] == 1e-2 + assert sb.scaling_factor[sb.flow_mol_phase["Liq"]] == 1e-2 + assert sb.scaling_factor[sb.flow_mol_phase["Vap"]] == 1e-2 + assert ( + sb.scaling_factor[ + sb.flow_mol_phase_comp["Liq", "benzene"] + ] + == 1e-2 + ) + assert ( + sb.scaling_factor[ + sb.flow_mol_phase_comp["Liq", "toluene"] + ] + == 1e-2 + ) + assert ( + sb.scaling_factor[ + sb.flow_mol_phase_comp["Vap", "benzene"] + ] + == 1e-2 + ) + assert ( + sb.scaling_factor[ + sb.flow_mol_phase_comp["Vap", "toluene"] + ] + == 1e-2 + ) + assert ( + sb.scaling_factor[sb.mole_frac_comp["benzene"]] + == 1000 + ) + assert ( + sb.scaling_factor[sb.mole_frac_comp["toluene"]] + == 1000 + ) + assert ( + sb.scaling_factor[ + sb.mole_frac_phase_comp["Liq", "benzene"] + ] + == 1000 + ) + assert ( + sb.scaling_factor[ + sb.mole_frac_phase_comp["Liq", "toluene"] + ] + == 1000 + ) + assert ( + sb.scaling_factor[ + sb.mole_frac_phase_comp["Vap", "benzene"] + ] + == 1000 + ) + assert ( + sb.scaling_factor[ + sb.mole_frac_phase_comp["Vap", "toluene"] + ] + == 1000 + ) + assert sb.scaling_factor[sb.pressure] == 1e-5 + assert sb.scaling_factor[sb.temperature] == 1e-2 + assert sb.scaling_factor[sb._teq["Vap", "Liq"]] == 1e-2 + assert sb.scaling_factor[sb._t1_Vap_Liq] == 1e-2 + + assert ( + sb.scaling_factor[ + sb._mole_frac_tbub["Vap", "Liq", "benzene"] + ] + == 1000 + ) + assert ( + sb.scaling_factor[ + sb._mole_frac_tbub["Vap", "Liq", "toluene"] + ] + == 1000 + ) + assert ( + sb.scaling_factor[ + sb._mole_frac_tdew["Vap", "Liq", "benzene"] + ] + == 1000 + ) + assert ( + sb.scaling_factor[ + sb._mole_frac_tdew["Vap", "Liq", "toluene"] + ] + == 1000 + ) + assert ( + sb.scaling_factor[sb.temperature_bubble["Vap", "Liq"]] + == 1e-2 + ) + assert ( + sb.scaling_factor[sb.temperature_dew["Vap", "Liq"]] + == 1e-2 + ) \ No newline at end of file diff --git a/property_packages/tests/test_compressor_pr.py b/property_packages/modular/tests/test_compressor_pr.py similarity index 74% rename from property_packages/tests/test_compressor_pr.py rename to property_packages/modular/tests/test_compressor_pr.py index 870a645..4fb7caf 100644 --- a/property_packages/tests/test_compressor_pr.py +++ b/property_packages/modular/tests/test_compressor_pr.py @@ -1,16 +1,15 @@ # Build and solve a heater block. -from ..build_package import build_package +from property_packages.build_package import build_package from pytest import approx # Import objects from pyomo package -from pyomo.environ import ConcreteModel, SolverFactory, value, units +from pyomo.environ import ConcreteModel, SolverFactory, value, units, Var # Import the main FlowsheetBlock from IDAES. The flowsheet block will contain the unit model from idaes.core import FlowsheetBlock from idaes.core.util.model_statistics import degrees_of_freedom from idaes.models.unit_models.heater import Heater from idaes.models.unit_models.pressure_changer import Pump -from idaes.core.util.tables import _get_state_from_port from idaes.models.unit_models.pressure_changer import PressureChanger, ThermodynamicAssumption def assert_approx(value, expected_value, error_margin): @@ -18,28 +17,28 @@ def assert_approx(value, expected_value, error_margin): tolerance = abs(percent_error * expected_value) assert approx(value, abs=tolerance) == expected_value + def test_compressor_asu(): m = ConcreteModel() m.fs = FlowsheetBlock(dynamic=False) m.fs.properties = build_package("peng-robinson", ["argon", "oxygen", "nitrogen"], ["Liq", "Vap"]) - + m.fs.compressor = PressureChanger( property_package=m.fs.properties, thermodynamic_assumption=ThermodynamicAssumption.isentropic, - compressor=True) + compressor=True + ) - m.fs.compressor.inlet.flow_mol[0].fix(1*units.kilomol/units.hour) - m.fs.compressor.inlet.pressure.fix(200000 * units.Pa) # doesn't work from 100,000 + m.fs.compressor.inlet.flow_mol[0].fix(1) + m.fs.compressor.inlet.pressure.fix(200000 * units.Pa) m.fs.compressor.inlet.temperature.fix((273.15+25) * units.K) - m.fs.compressor.inlet.mole_frac_comp[0, "argon"].fix(0.33) - m.fs.compressor.inlet.mole_frac_comp[0, "oxygen"].fix(0.33) - m.fs.compressor.inlet.mole_frac_comp[0, "nitrogen"].fix(0.33) + m.fs.compressor.inlet.mole_frac_comp[0, "argon"].fix(1/3) + m.fs.compressor.inlet.mole_frac_comp[0, "oxygen"].fix(1/3) + m.fs.compressor.inlet.mole_frac_comp[0, "nitrogen"].fix(1/3) - m.fs.compressor.outlet.pressure.fix(500000*units.Pa) + m.fs.compressor.outlet.pressure.fix(800000*units.Pa) m.fs.compressor.efficiency_isentropic.fix(0.75) - assert degrees_of_freedom(m) == 0 - m.fs.compressor.initialize() assert degrees_of_freedom(m) == 0 @@ -47,7 +46,8 @@ def test_compressor_asu(): solver = SolverFactory('ipopt') solver.solve(m, tee=True) - assert_approx(value(m.fs.compressor.outlet.temperature[0]), 466.12, 0.1) + assert_approx(value(m.fs.compressor.outlet.temperature[0]), 512.38, 0.1) + def test_expander_asu(): m = ConcreteModel() @@ -58,13 +58,13 @@ def test_expander_asu(): property_package=m.fs.properties, thermodynamic_assumption=ThermodynamicAssumption.adiabatic, compressor=False) - + m.fs.compressor.inlet.flow_mol[0].fix(1*units.kilomol/units.hour) m.fs.compressor.inlet.pressure.fix(1000000 * units.Pa) m.fs.compressor.inlet.temperature.fix((273.15+25) * units.K) - m.fs.compressor.inlet.mole_frac_comp[0, "argon"].fix(0.33) - m.fs.compressor.inlet.mole_frac_comp[0, "oxygen"].fix(0.33) - m.fs.compressor.inlet.mole_frac_comp[0, "nitrogen"].fix(0.33) + m.fs.compressor.inlet.mole_frac_comp[0, "argon"].fix(1/3) + m.fs.compressor.inlet.mole_frac_comp[0, "oxygen"].fix(1/3) + m.fs.compressor.inlet.mole_frac_comp[0, "nitrogen"].fix(1/3) m.fs.compressor.outlet.pressure.fix(100000*units.Pa) diff --git a/property_packages/tests/test_hc_pr.py b/property_packages/modular/tests/test_hc_pr.py similarity index 98% rename from property_packages/tests/test_hc_pr.py rename to property_packages/modular/tests/test_hc_pr.py index a9e5a4d..9944229 100644 --- a/property_packages/tests/test_hc_pr.py +++ b/property_packages/modular/tests/test_hc_pr.py @@ -38,11 +38,11 @@ ) from idaes.core.solvers import get_solver -from idaes.models.properties.modular_properties.state_definitions import FTPx +from idaes.models.properties.modular_properties.state_definitions import FTPx, FPhx from idaes.models.properties.modular_properties.phase_equil import SmoothVLE from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available -from ..build_package import build_package +from property_packages.build_package import build_package # ----------------------------------------------------------------------------- # Get default solver for testing @@ -163,7 +163,7 @@ def test_build(self): model.params.config.state_bounds, { "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), - "temperature": (13.95, 150, 350, pyunits.K), + "temperature": (1, 300, 3000, pyunits.K), "pressure": (5e4, 1e5, 1e6, pyunits.Pa), }, item_callback=_as_quantity, @@ -283,8 +283,8 @@ def test_build(self, model): assert isinstance(model.props[1].temperature, Var) assert value(model.props[1].temperature) == 295 - assert model.props[1].temperature.ub == 350 - assert model.props[1].temperature.lb == 13.95 + assert model.props[1].temperature.ub == 3000 + assert model.props[1].temperature.lb == 1 assert isinstance(model.props[1].mole_frac_comp, Var) assert len(model.props[1].mole_frac_comp) == 13 diff --git a/property_packages/tests/test_heat_exchanger_pr.py b/property_packages/modular/tests/test_heat_exchanger_pr.py similarity index 96% rename from property_packages/tests/test_heat_exchanger_pr.py rename to property_packages/modular/tests/test_heat_exchanger_pr.py index dc2281b..7f17345 100644 --- a/property_packages/tests/test_heat_exchanger_pr.py +++ b/property_packages/modular/tests/test_heat_exchanger_pr.py @@ -1,7 +1,6 @@ # Build and solve a heater block. -from ..build_package import build_package +from property_packages.build_package import build_package from pytest import approx - # Import objects from pyomo package from pyomo.environ import ConcreteModel, SolverFactory, value, units @@ -59,7 +58,7 @@ def test_heat_exchanger_bt(): solver = SolverFactory('ipopt') result = solver.solve(m) - assert value(m.fs.heat_exchanger.shell.properties_out[0].temperature) == approx(373.13, abs=1e-2) + assert value(m.fs.heat_exchanger.shell.properties_out[0].temperature) == approx(373.162, abs=1e-2) assert_approx(value(m.fs.heat_exchanger.tube.properties_out[0].temperature), 369.24, 0.5) def test_heat_exchanger_asu(): @@ -86,9 +85,9 @@ def test_heat_exchanger_asu(): assert degrees_of_freedom(m) == 8 m.fs.heat_exchanger.tube_inlet.flow_mol.fix(1*units.kilomol/units.hour) # mol/s - m.fs.heat_exchanger.tube_inlet.mole_frac_comp[0, "oxygen"].fix(0.33) - m.fs.heat_exchanger.tube_inlet.mole_frac_comp[0, "argon"].fix(0.33) - m.fs.heat_exchanger.tube_inlet.mole_frac_comp[0, "nitrogen"].fix(0.33) + m.fs.heat_exchanger.tube_inlet.mole_frac_comp[0, "oxygen"].fix(1/3) + m.fs.heat_exchanger.tube_inlet.mole_frac_comp[0, "argon"].fix(1/3) + m.fs.heat_exchanger.tube_inlet.mole_frac_comp[0, "nitrogen"].fix(1/3) m.fs.heat_exchanger.tube_inlet.pressure.fix(100000) # Pa m.fs.heat_exchanger.tube_inlet.temperature[0].fix((25 + 273.15)*units.K) @@ -104,4 +103,4 @@ def test_heat_exchanger_asu(): solver = SolverFactory('ipopt') result = solver.solve(m) - assert_approx(value(m.fs.heat_exchanger.shell.properties_out[0].temperature), 195.3 + 273.15, 0.1) \ No newline at end of file + assert_approx(value(m.fs.heat_exchanger.shell.properties_out[0].temperature), 195.3 + 273.15, 0.1) diff --git a/property_packages/tests/test_heater_pr.py b/property_packages/modular/tests/test_heater_pr.py similarity index 83% rename from property_packages/tests/test_heater_pr.py rename to property_packages/modular/tests/test_heater_pr.py index 06760e5..7cd4c3b 100644 --- a/property_packages/tests/test_heater_pr.py +++ b/property_packages/modular/tests/test_heater_pr.py @@ -1,7 +1,6 @@ # Build and solve a heater block. -from ..build_package import build_package +from property_packages.build_package import build_package from pytest import approx - # Import objects from pyomo package from pyomo.environ import ConcreteModel, SolverFactory, value, units @@ -9,7 +8,6 @@ from idaes.core import FlowsheetBlock from idaes.core.util.model_statistics import degrees_of_freedom from idaes.models.unit_models.heater import Heater -from idaes.core.util.tables import _get_state_from_port def assert_approx(value, expected_value, error_margin): percent_error = error_margin / 100 @@ -37,7 +35,7 @@ def test_heater_bt(): solver = SolverFactory('ipopt') solver.solve(m, tee=True) - assert_approx(value(_get_state_from_port(m.fs.heater.outlet,0).temperature), 363, 0.2) + assert_approx(value(m.fs.heater.outlet.temperature[0]), 363, 0.2) assert value(m.fs.heater.outlet.pressure[0]) == approx(101325) assert value(m.fs.heater.outlet.flow_mol[0]) == approx(1000/3600) @@ -49,9 +47,9 @@ def test_heater_asu(): m.fs.heater = Heater(property_package=m.fs.properties) m.fs.heater.inlet.flow_mol.fix(units.convert(1 * units.kilomol/units.hour, units.mol/units.s)) # mol/s - m.fs.heater.inlet.mole_frac_comp[0, "argon"].fix(0.33) - m.fs.heater.inlet.mole_frac_comp[0, "nitrogen"].fix(0.33) - m.fs.heater.inlet.mole_frac_comp[0, "oxygen"].fix(0.33) + m.fs.heater.inlet.mole_frac_comp[0, "argon"].fix(1/3) + m.fs.heater.inlet.mole_frac_comp[0, "nitrogen"].fix(1/3) + m.fs.heater.inlet.mole_frac_comp[0, "oxygen"].fix(1/3) m.fs.heater.inlet.pressure.fix(100000) m.fs.heater.inlet.temperature.fix(units.convert_temp_C_to_K(25)) m.fs.heater.outlet.temperature.fix(units.convert_temp_C_to_K(50)) @@ -68,9 +66,9 @@ def test_heater_asu(): m.fs.heater2 = Heater(property_package=m.fs.properties) m.fs.heater2.inlet.flow_mol.fix(units.convert(1 * units.kilomol/units.hour, units.mol/units.s)) # mol/s - m.fs.heater2.inlet.mole_frac_comp[0, "argon"].fix(0.33) - m.fs.heater2.inlet.mole_frac_comp[0, "nitrogen"].fix(0.33) - m.fs.heater2.inlet.mole_frac_comp[0, "oxygen"].fix(0.33) + m.fs.heater2.inlet.mole_frac_comp[0, "argon"].fix(1/3) + m.fs.heater2.inlet.mole_frac_comp[0, "nitrogen"].fix(1/3) + m.fs.heater2.inlet.mole_frac_comp[0, "oxygen"].fix(1/3) m.fs.heater2.inlet.pressure.fix(100000) m.fs.heater2.inlet.temperature.fix(units.convert_temp_C_to_K(25)) m.fs.heater2.outlet.temperature.fix(units.convert_temp_C_to_K(-15)) diff --git a/property_packages/tests/test_mixer_bt_pr.py b/property_packages/modular/tests/test_mixer_bt_pr.py similarity index 94% rename from property_packages/tests/test_mixer_bt_pr.py rename to property_packages/modular/tests/test_mixer_bt_pr.py index a0abc71..b6fc628 100644 --- a/property_packages/tests/test_mixer_bt_pr.py +++ b/property_packages/modular/tests/test_mixer_bt_pr.py @@ -1,5 +1,5 @@ # Build and solve a heater block. -from ..build_package import build_package +from property_packages.build_package import build_package from pytest import approx # Import objects from pyomo package @@ -9,7 +9,6 @@ from idaes.core import FlowsheetBlock from idaes.core.util.model_statistics import degrees_of_freedom from idaes.models.unit_models.heater import Heater -from idaes.core.util.tables import _get_state_from_port # Import objects from pyomo package from pyomo.environ import ConcreteModel, SolverFactory, value @@ -29,7 +28,6 @@ # Import the degrees_of_freedom function from the idaes.core.util.model_statistics package # DOF = Number of Model Variables - Number of Model Constraints from idaes.core.util.model_statistics import degrees_of_freedom -from idaes.core.util.tables import _get_state_from_port def assert_approx(value, expected_value, error_margin): percent_error = error_margin / 100 @@ -82,4 +80,4 @@ def test_mixer(): assert value(inlet_2.flow_mol_phase["Vap"]) == approx(0, abs=1e-2) assert value(inlet_2.flow_mol_phase["Liq"]) == approx(100, abs=1e-2) assert value(outlet.temperature) == approx(354.5, abs=1e-1) - assert value(outlet.pressure) == approx(101325*2, abs=1e-2) \ No newline at end of file + assert value(outlet.pressure) == approx(101325*2, abs=1e-2) diff --git a/property_packages/tests/test_pump_bt_pr.py b/property_packages/modular/tests/test_pump_bt_pr.py similarity index 88% rename from property_packages/tests/test_pump_bt_pr.py rename to property_packages/modular/tests/test_pump_bt_pr.py index 64e0604..a16948c 100644 --- a/property_packages/tests/test_pump_bt_pr.py +++ b/property_packages/modular/tests/test_pump_bt_pr.py @@ -1,7 +1,6 @@ # Build and solve a heater block. -from ..build_package import build_package +from property_packages.build_package import build_package from pytest import approx - # Import objects from pyomo package from pyomo.environ import ConcreteModel, SolverFactory, value, units @@ -10,7 +9,6 @@ from idaes.core.util.model_statistics import degrees_of_freedom from idaes.models.unit_models.heater import Heater from idaes.models.unit_models.pressure_changer import Pump -from idaes.core.util.tables import _get_state_from_port def assert_approx(value, expected_value, error_margin): percent_error = error_margin / 100 @@ -40,6 +38,6 @@ def test_pump(): solver = SolverFactory('ipopt') solver.solve(m, tee=True) - assert_approx(value(_get_state_from_port(m.fs.pump.outlet,0).temperature), 353, 0.2) + assert_approx(value(m.fs.pump.outlet.temperature[0]), 353, 0.2) assert value(m.fs.pump.outlet.pressure[0]) == approx(201325) assert value(m.fs.pump.outlet.flow_mol[0]) == approx(100) \ No newline at end of file diff --git a/property_packages/modular/tests/test_sb_constraints.py b/property_packages/modular/tests/test_sb_constraints.py new file mode 100644 index 0000000..ab2888d --- /dev/null +++ b/property_packages/modular/tests/test_sb_constraints.py @@ -0,0 +1,196 @@ +""" +Test initializing and solving a state block with constraints +rather than fixing the state variables directly. +""" + +from pytest import approx +from pyomo.environ import ( + ConcreteModel, + SolverFactory, + value, + units, + assert_optimal_termination, +) +from pyomo.core.base.constraint import Constraint, ScalarConstraint +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import degrees_of_freedom +from ..template_builder import build_config + + +def flowsheet(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = build_config( + "peng-robinson", ["benzene", "toluene"], ["Liq", "Vap"] + ) + return m + + +def build_state(m): + m.fs.state = m.fs.properties.build_state_block([0], defined_state=True) + return m.fs.state[0] + + +def initialize(m): + assert degrees_of_freedom(m) == 0 + m.fs.state.initialize() + assert degrees_of_freedom(m) == 0 + + +def solve(m): + solver = SolverFactory("ipopt") + res = solver.solve(m, tee=True) + assert_optimal_termination(res) + + +def test_constrain(): + m = flowsheet() + sb = build_state(m) + + # fix flow_mol by string + c = sb.constrain("flow_mol", 10) + assert c == sb.flow_mol + assert c.value == 10 + assert c.is_fixed() + + # constrain flow_mass by string + c = sb.constrain("flow_mass", 10) + assert type(c) == ScalarConstraint + assert c in sb.component_data_objects(Constraint) + assert getattr(sb.constraints, "flow_mass") == c + + +def test_constrain_component(): + m = flowsheet() + sb = build_state(m) + + # fix flow_mol by component + c = sb.constrain_component(sb.flow_mol, 10) + assert c == sb.flow_mol + assert c.value == 10 + assert c.is_fixed() + + # constrain flow_mass by component + c = sb.constrain_component(sb.flow_mass, 10) + assert type(c) == ScalarConstraint + assert c in sb.component_data_objects(Constraint) + assert getattr(sb.constraints, "flow_mass") == c + + +def test_state_vars(): + # default state vars: flow_mol, temperature, pressure, mole_frac_comp + m = flowsheet() + sb = build_state(m) + sb.constrain("flow_mol", 1) + sb.constrain("temperature", 290) + sb.constrain("pressure", 100000) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + initialize(m) + solve(m) + assert value(sb.temperature) == approx(290) + assert value(sb.pressure) == approx(100000) + assert value(sb.flow_mol) == approx(1) + assert value(sb.mole_frac_comp["benzene"]) == approx(0.5) + assert value(sb.mole_frac_comp["toluene"]) == approx(0.5) + + +def test_enth_mol(): + m = flowsheet() + sb = build_state(m) + sb.constrain("flow_mol", 1) + sb.constrain("enth_mol", 30611.284116732746) # J/mol + sb.constrain("pressure", 100000) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + initialize(m) + solve(m) + assert value(sb.temperature) == approx(290) + assert value(sb.enth_mol) == approx(30611.284116732746) + + +def test_enth_mass(): + m = flowsheet() + sb = build_state(m) + sb.constrain("flow_mass", 1) + sb.constrain("enth_mass", 359603.35471696255) + sb.constrain("pressure", 101325) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + initialize(m) + solve(m) + assert value(sb.temperature) == approx(290, rel=1e-3) + assert value(sb.enth_mass) == approx(359603.35471696255) + + +def test_entr_mol(): + m = flowsheet() + sb = build_state(m) + sb.constrain("flow_mass", 1) + sb.constrain("entr_mol", -388.1271856851202) + sb.constrain("pressure", 101325) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + initialize(m) + solve(m) + assert value(sb.temperature) == approx(290) + assert value(sb.entr_mol) == approx(-388.1271856851202) + + +def test_entr_mass(): + m = flowsheet() + sb = build_state(m) + sb.constrain("flow_mass", 1) + sb.constrain("entr_mass", -4559.489810913312) + sb.constrain("pressure", 101325) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + initialize(m) + solve(m) + assert value(sb.temperature) == approx(290) + assert value(sb.entr_mass) == approx(-4559.489810913312) + + +def test_flow_mass(): + m = flowsheet() + sb = build_state(m) + sb.constrain("flow_mass", 0.08512513500000002) + sb.constrain("temperature", 290) + sb.constrain("pressure", 101325) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + initialize(m) + solve(m) + assert value(sb.flow_mass) == approx(0.08512513500000002) + assert value(sb.flow_mol) == approx(1) + + +def test_flow_vol(): + m = flowsheet() + sb = build_state(m) + sb.constrain("flow_vol", 9.660626354419042e-05) + sb.constrain("temperature", 290) + sb.constrain("pressure", 101325) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + initialize(m) + solve(m) + assert value(sb.flow_vol) == approx(9.660626354419042e-05) + assert value(sb.flow_mol) == approx(1, rel=1e-3) + + +def test_enthalpy_temperature(): + # use enthalpy and temperature to define pressure + # TODO: this test is not currently working + # related: https://github.com/IDAES/idaes-pse/pull/1554 + m = flowsheet() + sb = build_state(m) + sb.constrain("flow_mass", 1) + sb.constrain("temperature", 290) + sb.constrain("enth_mol", 30611.284116732746) + sb.mole_frac_comp["benzene"].fix(0.5) + sb.mole_frac_comp["toluene"].fix(0.5) + return # fails to initialize + initialize(m) + solve(m) + assert value(sb.pressure) == approx(100000) diff --git a/property_packages/tests/test_separator_H20_C02_pr.py b/property_packages/modular/tests/test_separator_H20_C02_pr.py similarity index 97% rename from property_packages/tests/test_separator_H20_C02_pr.py rename to property_packages/modular/tests/test_separator_H20_C02_pr.py index 415c0c9..476b28c 100644 --- a/property_packages/tests/test_separator_H20_C02_pr.py +++ b/property_packages/modular/tests/test_separator_H20_C02_pr.py @@ -31,10 +31,8 @@ from pyomo.environ import ConcreteModel, SolverFactory, value, units -from idaes.core.util.tables import _get_state_from_port - # Import the build function to create a property package -from ..build_package import build_package +from property_packages.build_package import build_package from pytest import approx diff --git a/property_packages/modular/tests/test_state_bounds.py b/property_packages/modular/tests/test_state_bounds.py new file mode 100644 index 0000000..a17a786 --- /dev/null +++ b/property_packages/modular/tests/test_state_bounds.py @@ -0,0 +1,76 @@ +""" +Try to test the property package at the bounds +of the state vars. For example, going to the upper +bound of the temperature or pressure. +""" + +from pyomo.environ import ConcreteModel, SolverFactory, check_optimal_termination +from idaes.core import FlowsheetBlock +from idaes.models.unit_models.pressure_changer import PressureChanger +from idaes.models.unit_models.heater import Heater +from property_packages.build_package import build_package + + +def test_compressor(): + """ + Test a compressor with pressure change from lb to ub + """ + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = build_package( + "peng-robinson", + ["methane", "ethane", "propane", "carbon dioxide", "nitrogen"], + ["Liq", "Vap"], + ) + m.fs.compressor = PressureChanger( + property_package=m.fs.properties, + thermodynamic_assumption="isentropic", + compressor=True, + ) + m.fs.compressor.efficiency_isentropic.fix(0.75) + m.fs.compressor.inlet.flow_mol.fix(1) + m.fs.compressor.inlet.pressure.fix(m.fs.compressor.inlet.pressure[0].lb) + m.fs.compressor.inlet.temperature.fix(300) + m.fs.compressor.inlet.mole_frac_comp[0, "methane"].fix(0.87) + m.fs.compressor.inlet.mole_frac_comp[0, "ethane"].fix(0.06) + m.fs.compressor.inlet.mole_frac_comp[0, "propane"].fix(0.03) + m.fs.compressor.inlet.mole_frac_comp[0, "carbon dioxide"].fix(0.02) + m.fs.compressor.inlet.mole_frac_comp[0, "nitrogen"].fix(0.02) + m.fs.compressor.outlet.pressure.fix(m.fs.compressor.outlet.pressure[0].ub) + + m.fs.compressor.initialize(outlvl=1) + opt = SolverFactory("ipopt") + res = opt.solve(m, tee=True) + assert check_optimal_termination(res) + + +def test_heater(): + """ + Test a heater with temperature change from lb to ub + """ + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = build_package( + "peng-robinson", + ["nitrogen", "carbon dioxide", "water"], + ["Liq", "Vap"], + ) + m.fs.heater = Heater( + property_package=m.fs.properties, + has_pressure_change=False, + ) + m.fs.heater.inlet.flow_mol.fix(1) + m.fs.heater.inlet.pressure.fix(100000) + # temperature currently doesn't have a good method for determining lb + # we will do 150 K for now + m.fs.heater.inlet.temperature.fix(150) + m.fs.heater.inlet.mole_frac_comp[0, "nitrogen"].fix(0.79) + m.fs.heater.inlet.mole_frac_comp[0, "carbon dioxide"].fix(0.2) + m.fs.heater.inlet.mole_frac_comp[0, "water"].fix(0.01) + m.fs.heater.outlet.temperature.fix(m.fs.heater.outlet.temperature[0].ub) + + m.fs.heater.initialize(outlvl=1) + opt = SolverFactory("ipopt") + res = opt.solve(m, tee=True) + + assert check_optimal_termination(res) diff --git a/property_packages/modular/tests/test_state_definitions.py b/property_packages/modular/tests/test_state_definitions.py deleted file mode 100644 index eb64b30..0000000 --- a/property_packages/modular/tests/test_state_definitions.py +++ /dev/null @@ -1,65 +0,0 @@ -from ..template_builder import build_config -from pytest import approx -from pyomo.environ import ConcreteModel, SolverFactory, value, units -from idaes.core import FlowsheetBlock -from idaes.models.unit_models import Compressor - -def flowsheet(): - m = ConcreteModel() - m.fs = FlowsheetBlock(dynamic=False) - m.fs.properties = build_config("peng-robinson",["benzene","toluene"],["Liq","Vap"]) - return m - -def solve(m): - solver = SolverFactory('ipopt') - solver.solve(m, tee=True) - -def test_state_definition2(): - m = flowsheet() - m.fs.sb = m.fs.properties.build_state_block(m.fs.time) - sb = m.fs.sb[0] - sb.constrain("temperature", 280) - sb.constrain("pressure", 101325) - sb.constrain("flow_mass", 1) - sb.mole_frac_comp["benzene"].fix(0.5) - sb.mole_frac_comp["toluene"].fix(0.5) - m.fs.sb.initialize() - solve(m) - assert value(sb.temperature) == approx(280) - - -def test_state_definition_temp(): - m = flowsheet() - m.fs.sb = m.fs.properties.build_state_block(m.fs.time) - sb = m.fs.sb[0] - sb.constrain("smooth_temperature", 273.5809) - sb.constrain("pressure", 101325) - sb.constrain("flow_mass", 1) - sb.mole_frac_comp["benzene"].fix(0.5) - sb.mole_frac_comp["toluene"].fix(0.5) - m.fs.sb.initialize() - solve(m) - assert value(sb.enth_mass) == approx(1878.712) - - -def test_initialise_compressor(): - # The purpose of this test is to make sure the compressor can initialise. - # We have had some errors with Too Few Degrees of Freedom - # being thrown in initialisation. This is because constraints are - # being fixed in the inlet state block instead of variables. - # This asserts that the appropriate variables are unfixed. - m = flowsheet() - m.fs.compressor = Compressor(property_package=m.fs.properties) - inlet = m.fs.compressor.control_volume.properties_in[0] - outlet = m.fs.compressor.control_volume.properties_out[0] - inlet.constrain("smooth_temperature", 273.5809) - inlet.constrain("pressure", 101325) - inlet.constrain("flow_mass", 1) - inlet.mole_frac_comp["benzene"].fix(0.5) - inlet.mole_frac_comp["toluene"].fix(0.5) - m.fs.compressor.deltaP.fix(100000) - m.fs.compressor.initialize() - solve(m) - assert value(outlet.temperature) == approx(393.5689573) - - diff --git a/property_packages/tests/test_bt_pr.py b/property_packages/tests/test_bt_pr.py deleted file mode 100644 index a6ba472..0000000 --- a/property_packages/tests/test_bt_pr.py +++ /dev/null @@ -1,438 +0,0 @@ -# Build and solve a state block. -from ..build_package import build_package -from pytest import approx - -# Import objects from pyomo package -from pyomo.environ import ConcreteModel, value - -# Import the main FlowsheetBlock from IDAES. The flowsheet block will contain the unit model -from idaes.core import FlowsheetBlock -from idaes.core.solvers import get_solver - -from pyomo.util.check_units import assert_units_consistent - -from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available -from pyomo.environ import check_optimal_termination, ConcreteModel, Objective -import idaes.core.util.scaling as iscale - -solver = get_solver(solver="ipopt") - -import idaes.logger as idaeslog -SOUT = idaeslog.INFO - -""" -Test Suite Sourced From IDAES - -URL: https://github.com/IDAES/idaes-pse/blob/41bb3c9728ea227fa8bb0fa1bf35b8deec467783/idaes/models/ -properties/modular_properties/examples/tests/test_BT_PR.py -""" - -""" ---------------------------------------------- -Goal for absolute error margins ---------------------------------------------- -- temp, pressure, specific volume : 0.5% -- compositional things : 2% (mole frac, flow, etc) -- thermodynamic values: 5% (entr, enth) ---------------------------------------------- -""" - -""" -Helper function to calculate the absolute error margin of a -value and asserts whether or not value lies within range -- Accepts: percent_error (as whole percent) -- Returns: None -""" -def assert_approx(value, expected_value, error_margin): - percent_error = error_margin / 100 - tolerance = abs(percent_error * expected_value) - assert approx(value, abs=tolerance) == expected_value - -def get_m(): - m = ConcreteModel() - m.fs = FlowsheetBlock(dynamic=False) - m.fs.props = build_package("peng-robinson", ["benzene", "toluene"], ["Liq", "Vap"]) - m.fs.state = m.fs.props.build_state_block([1], defined_state=True) - iscale.calculate_scaling_factors(m.fs.props) - iscale.calculate_scaling_factors(m.fs.state[1]) - return m - -def test_T_sweep(): - m = get_m() - - assert_units_consistent(m) - - m.fs.obj = Objective(expr=(m.fs.state[1].temperature - 510) ** 2) - m.fs.state[1].temperature.setub(600) - - for logP in [9.5, 10, 10.5, 11, 11.5, 12]: - m.fs.obj.deactivate() - - m.fs.state[1].flow_mol.fix(100) - m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) - m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) - m.fs.state[1].temperature.fix(300) - m.fs.state[1].pressure.fix(10 ** (0.5 * logP)) - - m.fs.state.initialize() - - m.fs.state[1].temperature.unfix() - m.fs.obj.activate() - - results = solver.solve(m) - - assert check_optimal_termination(results) - assert m.fs.state[1].flow_mol_phase["Liq"].value <= 1e-2 - -def test_P_sweep(): - m = get_m() - - for T in range(370, 500, 25): - m.fs.state[1].flow_mol.fix(100) - m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) - m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) - m.fs.state[1].temperature.fix(T) - m.fs.state[1].pressure.fix(1e5) - print(T) - m.fs.state.initialize() - - results = solver.solve(m) - - assert check_optimal_termination(results) - - while m.fs.state[1].pressure.value <= 1e6: - - results = solver.solve(m) - assert check_optimal_termination(results) - - m.fs.state[1].pressure.value = m.fs.state[1].pressure.value + 1e5 - -def test_T350_P1_x5(): - m = get_m() - - m.fs.state[1].flow_mol.fix(100) - m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) - m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) - m.fs.state[1].temperature.fix(350) - m.fs.state[1].pressure.fix(1e5) - - # Trigger build of enthalpy and entropy - m.fs.state[1].enth_mol_phase - m.fs.state[1].entr_mol_phase - - m.fs.state.initialize(outlvl=SOUT) - - results = solver.solve(m) - - # Check for optimal solution - assert check_optimal_termination(results) - - assert_approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 365, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Liq"]), 0.0035346, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Vap"]), 0.966749, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 0.894676, 2) # Investigate - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 0.347566, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 0.971072, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 0.959791, 0.5) - - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 0.5, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 0.5, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 0.70584, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 0.29416, 2) - - assert_approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 38942.8, 5) - assert_approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 78048.7, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Liq"]), -361.794, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Vap"]), -264.0181, 5) - -def test_T350_P5_x5(): - m = get_m() - - m.fs.state[1].flow_mol.fix(100) - m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) - m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) - m.fs.state[1].temperature.fix(350) - m.fs.state[1].pressure.fix(5e5) - - # Trigger build of enthalpy and entropy - m.fs.state[1].enth_mol_phase - m.fs.state[1].entr_mol_phase - - m.fs.state.initialize(outlvl=SOUT) - - results = solver.solve(m) - - # Check for optimal solution - assert check_optimal_termination(results) - - assert_approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 431.47, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Liq"]), 0.01766, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Vap"]), 0.80245, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 0.181229, 1.5) # Investigate - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 0.070601, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 0.856523, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 0.799237, 0.5) - - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 0.5, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 0.5, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 0.65415, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 0.34585, 2) - - assert_approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 38966.9, 5) - assert_approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 75150.7, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Liq"]), -361.8433, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Vap"]), -281.9703, 5) - -def test_T450_P1_x5(): - m = get_m() - - m.fs.state[1].flow_mol.fix(100) - m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) - m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) - m.fs.state[1].temperature.fix(450) - m.fs.state[1].pressure.fix(1e5) - - # Trigger build of enthalpy and entropy - m.fs.state[1].enth_mol_phase - m.fs.state[1].entr_mol_phase - - m.fs.state.initialize(outlvl=SOUT) - - results = solver.solve(m) - - # Check for optimal solution - assert check_optimal_termination(results) - - assert_approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 371.4, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Liq"]), 0.0033583, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Vap"]), 0.9821368, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 8.069323, 1) # Investigate - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 4.304955, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 0.985365, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 0.979457, 0.5) - - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 0.29861, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 0.70139, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 0.5, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 0.5, 2) - - assert_approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 49441.2, 5) - assert_approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 84175.1, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Liq"]), -328.766, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Vap"]), -241.622, 5) - -def test_T450_P5_x5(): - m = get_m() - - m.fs.state[1].flow_mol.fix(100) - m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) - m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) - m.fs.state[1].temperature.fix(450) - m.fs.state[1].pressure.fix(5e5) - - # Trigger build of enthalpy and entropy - m.fs.state[1].enth_mol_phase - m.fs.state[1].entr_mol_phase - - m.fs.state.initialize(outlvl=SOUT) - - results = solver.solve(m) - - # Check for optimal solution - assert check_optimal_termination(results) - - assert_approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 436.93, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Liq"]), 0.0166181, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Vap"]), 0.9053766, 0.5) - - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1.63308, 1) # investigate - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 0.873213, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 0.927534, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 0.898324, 0.5) - - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 0.3488737, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 0.6511263, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 0.5, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 0.5, 2) - - assert_approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 51095.2, 5) - assert_approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 83362.3, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Liq"]), -326.299, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Vap"]), -256.198, 5) - -def test_T368_P1_x5(): - m = get_m() - - m.fs.state[1].flow_mol.fix(100) - m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) - m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) - m.fs.state[1].temperature.fix(368) - m.fs.state[1].pressure.fix(1e5) - - # Trigger build of enthalpy and entropy - m.fs.state[1].enth_mol_phase - m.fs.state[1].entr_mol_phase - - m.fs.state.initialize(outlvl=SOUT) - - results = solver.solve(m) - - # Check for optimal solution - assert check_optimal_termination(results) - - assert_approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 368, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Liq"]), 0.003504, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Vap"]), 0.97, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1.492049, 1.5) # Investigate - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 0.621563, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 0.97469, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 0.964642, 0.5) - - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 0.4012128, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 0.5987872, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 0.6141738, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 0.3858262, 2) - - assert_approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 38235.1, 5) - assert_approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 77155.4, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Liq"]), -359.256, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Vap"]), -262.348, 5) - -def test_T376_P1_x2(): - m = get_m() - - m.fs.state[1].flow_mol.fix(100) - m.fs.state[1].mole_frac_comp["benzene"].fix(0.2) - m.fs.state[1].mole_frac_comp["toluene"].fix(0.8) - m.fs.state[1].temperature.fix(376) - m.fs.state[1].pressure.fix(1e5) - - # Trigger build of enthalpy and entropy - m.fs.state[1].enth_mol_phase - m.fs.state[1].entr_mol_phase - - m.fs.state.initialize(outlvl=SOUT) - - results = solver.solve(m) - - # Check for optimal solution - assert check_optimal_termination(results) - - assert_approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 376, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Liq"]), 0.00361333, 0.5) - assert_approx(value(m.fs.state[1].compress_fact_phase["Vap"]), 0.968749, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1.8394188, 1.5) # Investigate - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 0.7871415, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 0.9763608, 0.5) - assert_approx(value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 0.9663611, 0.5) - - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 0.17342, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 0.82658, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 0.3267155, 2) - assert_approx(value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 0.6732845, 2) - - assert_approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 31535.8, 5) - assert_approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 69175.3, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Liq"]), -369.033, 5) - assert_approx(value(m.fs.state[1].entr_mol_phase["Vap"]), -273.513, 5) - -def test_basic_scaling(): - m = get_m() - - assert len(m.fs.state[1].scaling_factor) == 23 - assert m.fs.state[1].scaling_factor[m.fs.state[1].flow_mol] == 1e-2 - assert m.fs.state[1].scaling_factor[m.fs.state[1].flow_mol_phase["Liq"]] == 1e-2 - assert m.fs.state[1].scaling_factor[m.fs.state[1].flow_mol_phase["Vap"]] == 1e-2 - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].flow_mol_phase_comp["Liq", "benzene"] - ] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].flow_mol_phase_comp["Liq", "toluene"] - ] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].flow_mol_phase_comp["Vap", "benzene"] - ] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].flow_mol_phase_comp["Vap", "toluene"] - ] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[m.fs.state[1].mole_frac_comp["benzene"]] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[m.fs.state[1].mole_frac_comp["toluene"]] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"] - ] - == 1000 - ) - assert m.fs.state[1].scaling_factor[m.fs.state[1].pressure] == 1e-5 - assert m.fs.state[1].scaling_factor[m.fs.state[1].temperature] == 1e-2 - assert m.fs.state[1].scaling_factor[m.fs.state[1]._teq["Vap", "Liq"]] == 1e-2 - assert m.fs.state[1].scaling_factor[m.fs.state[1]._t1_Vap_Liq] == 1e-2 - - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1]._mole_frac_tbub["Vap", "Liq", "benzene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1]._mole_frac_tbub["Vap", "Liq", "toluene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1]._mole_frac_tdew["Vap", "Liq", "benzene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1]._mole_frac_tdew["Vap", "Liq", "toluene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[m.fs.state[1].temperature_bubble["Vap", "Liq"]] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[m.fs.state[1].temperature_dew["Vap", "Liq"]] - == 1e-2 - ) \ No newline at end of file diff --git a/property_packages/tests/test_build.py b/property_packages/tests/test_build.py deleted file mode 100644 index 1fd4dfd..0000000 --- a/property_packages/tests/test_build.py +++ /dev/null @@ -1,33 +0,0 @@ -import pytest -from pyomo.environ import ( - assert_optimal_termination, - ConcreteModel, - Set, - value, - Var, - units as pyunits, -) -from pyomo.util.check_units import assert_units_consistent -import pyomo.common.unittest as unittest - -from idaes.core import Component -from idaes.core.util.model_statistics import ( - degrees_of_freedom, - fixed_variables_set, - activated_constraints_set, -) -from idaes.core.solvers import get_solver - -from idaes.models.properties.modular_properties.state_definitions import FTPx -from idaes.models.properties.modular_properties.phase_equil import SmoothVLE - -from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available -from ..build_package import build_package - -def test_build(): - build_package("helmholtz", ["water"]) - build_package("peng-robinson", ["benzene", "toluene"]) - build_package("helmholtz", ["water"], ["Liq"]) - build_package("peng-robinson", ["benzene", "toluene"], ["Liq"]) - build_package("helmholtz", ["water"], ["Vap", "Liq"]) - build_package("peng-robinson", ["benzene", "toluene"], ["Vap", "Liq"]) \ No newline at end of file diff --git a/property_packages/utils/add_extra_expressions.py b/property_packages/utils/add_extra_expressions.py index effb40b..c726203 100644 --- a/property_packages/utils/add_extra_expressions.py +++ b/property_packages/utils/add_extra_expressions.py @@ -4,10 +4,6 @@ def add_extra_expressions(sb: Block): """ idaes state blocks don't support all the properties we need, so we add some extra expressions here """ - #if not hasattr(sb, "smooth_pressure"): - # sb.add_component("smooth_pressure", Expression(expr=sb.pressure + value(sb.enth_mol)*0.0001 * units.Pa)) - if not hasattr(sb, "smooth_temperature"): - sb.add_component("smooth_temperature", Expression(expr=sb.temperature + (sb.enth_mol / (1 * units.J/units.mol))*0.000001 * units.K)) if not hasattr(sb, "enth_mass"): sb.add_component("enth_mass", Expression(expr=(sb.flow_mol * sb.enth_mol) / sb.flow_mass)) if not hasattr(sb, "entr_mass"): diff --git a/property_packages/utils/fix_state_vars.py b/property_packages/utils/fix_state_vars.py new file mode 100644 index 0000000..d05c9e8 --- /dev/null +++ b/property_packages/utils/fix_state_vars.py @@ -0,0 +1,96 @@ +from idaes.core.util.exceptions import ConfigurationError + +""" +This method is analogous to the fix_state_vars method in +idaes.core.util.initialization, but it allows us to handle +the case where we have constraints that define the state. +""" + +def fix_state_vars(blk, state_args=None): + """ + Method for fixing state variables within StateBlocks. Method takes an + optional argument of values to use when fixing variables. + + Args: + blk : An IDAES StateBlock object in which to fix the state variables + state_args : a dict containing values to use when fixing state + variables. Keys must match with names used in the + define_state_vars method, and indices of any variables must + agree. + + Returns: + A dict keyed by block index, state variable name (as defined by + define_state_variables) and variable index indicating the fixed status + of each variable before the fix_state_vars method was applied. + """ + if state_args is None: + state_args = {} + + flags = {} + for k, b in blk.items(): + flow_exprs = { + "flow_mol", + "flow_mass", + "flow_vol" + } + other_exprs = { + "pressure", + "enth_mol", + "enth_mass", + "entr_mol", + "entr_mass", + "temperature", + "total_energy_flow", + "custom_vapor_frac", + } + for n, v in b.define_state_vars().items(): + fix_var = True + if n in flow_exprs: + for prop_name in flow_exprs: + if hasattr(b.constraints, prop_name): + fix_var = False + break + elif n in other_exprs: + if not v.is_fixed(): + # check if any of the constraints exist + for prop_name in other_exprs: + if hasattr(b.constraints, prop_name): + # don't fix this variable - it is defined by a constraint + other_exprs.remove(prop_name) + fix_var = False + break + for i in v: + flags[k, n, i] = v[i].is_fixed() + if fix_var and not v[i].is_fixed(): + # If not fixed, fix at either guess provided or current value + if n in state_args: + # Try to get initial guess from state_args + try: + if i is None: + val = state_args[n] + else: + val = state_args[n][i] + except KeyError: + raise ConfigurationError( + "Indexes in state_args did not agree with " + "those of state variable {}. Please ensure " + "that indexes for initial guesses are correct.".format( + n + ) + ) + v[i].fix(val) + else: + # No guess, try to use current value + if v[i].value is not None: + v[i].fix() + else: + # No initial value - raise Exception before this + # gets to a solver. + raise ConfigurationError( + "State variable {} does not have a value " + "assigned. This usually occurs when a Var " + "is not assigned an initial value when it is " + "created. Please ensure all variables have " + "valid values before fixing them.".format(v.name) + ) + return flags \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 668af92..24e8ab2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ build-backend = "hatchling.build" [project] name = "ahuora-compounds" -version = "0.0.5" +version = "0.0.10" authors = [ { name="Example Author", email="author@example.com" }, @@ -34,4 +34,4 @@ property_packages = [ Homepage = "https://ahuora.org.nz" [tool.hatch.build.targets.wheel] -packages = ["compounds","property_packages"] \ No newline at end of file +packages = ["compounds","property_packages"] diff --git a/test.py b/test.py deleted file mode 100644 index 6b9ed89..0000000 --- a/test.py +++ /dev/null @@ -1,88 +0,0 @@ -### Imports -from pyomo.environ import ConcreteModel, SolverFactory, SolverStatus, TerminationCondition, Block, TransformationFactory, Constraint, value -from pyomo.network import SequentialDecomposition, Port, Arc -from pyomo.core.base.units_container import _PyomoUnit, units as pyomo_units -from idaes.core import FlowsheetBlock -from idaes.core.util.model_statistics import report_statistics, degrees_of_freedom, fixed_variables_set, activated_block_component_generator -import idaes.logger as idaeslog -from idaes.models.properties.general_helmholtz import HelmholtzParameterBlock, PhaseType, StateVars, AmountBasis -from idaes.models.properties.iapws95 import Iapws95ParameterBlock -from idaes.models.unit_models.heater import Heater - - -### Utility Methods -def units(item: str) -> _PyomoUnit: - ureg = pyomo_units._pint_registry - pint_unit = getattr(ureg, item) - return _PyomoUnit(pint_unit, ureg) - - -### Build Model -m = ConcreteModel() -m.fs = FlowsheetBlock(dynamic=False) - -# Set up property packages -m.fs.properties = HelmholtzParameterBlock( - pure_component="h2o", - phase_presentation=PhaseType.MIX, - state_vars=StateVars.PH, - amount_basis=AmountBasis.MOLE, -) -# m.fs.properties = Iapws95ParameterBlock() -# Create unit models - -# HT-772 -m.fs.HT_772 = Heater( - property_package=m.fs.properties, - has_pressure_change=True, -) -inlet_state = m.fs.HT_772.control_volume.properties_in[0] -outlet_state = m.fs.HT_772.control_volume.properties_out[0] -inlet_state.flow_mol.fix(1.0 * units("mol/s")) -inlet_state.pressure.fix(100.0 * units("kPa")) -# inlet_state.temperature.fix(283.15 * units("K")) -# inlet_state.vapor_frac.fix(0.0) -# inlet_state.enth_mol.fix(80000 * units("J/mol")) -# m.fs.HT_772.deltaP[0].fix(-0.00001 * (inlet_state.enth_mol - outlet_state.enth_mol)) -# outlet_state.pressure.fix(99.999 * units("kPa")) -# outlet_state.enth_mol.fix(80000 * units("J/mol")) -# outlet_state.temperature.fix(293.15 * units("K")) -# outlet_state.enth_mol.fix(1513.4 * units("J/mol")) -m.constraint0 = Constraint(expr=inlet_state.temperature == 383.15) # K -m.contraint1 = Constraint(expr=outlet_state.temperature == 393.15) # K -m.fs.HT_772.constraint2 = Constraint(expr=m.fs.HT_772.deltaP[0] == (0.0001 * (outlet_state.enth_mol - inlet_state.enth_mol))) -print(hasattr(m.fs.HT_772, "deltaT")) - - -### Check Model Status -report_statistics(m) -print("Degrees of freedom:", degrees_of_freedom(m)) - - -### Initialize Model -m.fs.HT_772.initialize(outlvl=idaeslog.INFO) - -outlet_state.enth_mol.unfix() - -m.fs.HT_772.report() - -fixed = fixed_variables_set(m) -for v in fixed: - print(v, "fixed to", value(v)) -constraints = activated_block_component_generator(m, Constraint) -for c in constraints: - print(c, "expression:", c.expr) - - -### Solve -solver = SolverFactory("ipopt") -result = solver.solve(m, tee=True) -# check for optimal solution -if result.solver.status != SolverStatus.ok or result.solver.termination_condition != TerminationCondition.optimal: - raise Exception("Solver did not converge to an optimal solution") - - -### Report -m.fs.HT_772.report() - -#print(value(m.fs.HT_772.deltaP[0])) \ No newline at end of file