diff --git a/idaes/core/base/control_volume0d.py b/idaes/core/base/control_volume0d.py index 9ece9ef157..4cc05b1f22 100644 --- a/idaes/core/base/control_volume0d.py +++ b/idaes/core/base/control_volume0d.py @@ -68,6 +68,8 @@ class ControlVolume0DScaler(ControlVolumeScalerBase): "phase_fraction": 10, # May have already been created by property package } + _state_block_ref = "properties_out" + def variable_scaling_routine( self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None ): diff --git a/idaes/core/base/control_volume1d.py b/idaes/core/base/control_volume1d.py index 03889f2d34..985123fc9d 100644 --- a/idaes/core/base/control_volume1d.py +++ b/idaes/core/base/control_volume1d.py @@ -25,6 +25,7 @@ # Import Pyomo libraries from pyomo.environ import ( + ComponentMap, Constraint, Expression, Param, @@ -47,8 +48,11 @@ MaterialFlowBasis, MaterialBalanceType, ) +from idaes.core.base.control_volume_base import ControlVolumeScalerBase +from idaes.core.scaling import DefaultScalingRecommendation from idaes.core.util.exceptions import ( BalanceTypeNotSupportedError, + BurntToast, ConfigurationError, PropertyNotSupportedError, ) @@ -58,7 +62,7 @@ import idaes.logger as idaeslog -__author__ = "Andrew Lee, Jaffer Ghouse" +__author__ = "Andrew Lee, Jaffer Ghouse, Douglas Allan" _log = idaeslog.getLogger(__name__) @@ -77,6 +81,258 @@ class DistributedVars(Enum): uniform = 1 +class ControlVolume1DScaler(ControlVolumeScalerBase): + """ + Scaler object for the ControlVolume1D + """ + + DEFAULT_SCALING_FACTORS = { + # We could scale length and area by magnitude if it were + # being fixed by the user, but we often have the volume + # given by an equality constraint involving geometry in + # the parent unit model. + "area": DefaultScalingRecommendation.userInputRequired, + "length": DefaultScalingRecommendation.userInputRequired, + "phase_fraction": 10, # May have already been created by property package + } + + _state_block_ref = "properties" + _weight_attr_name = "length" + + def variable_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None + ): + """ + Routine to apply scaling factors to variables in model. + + Derived classes must overload this method. + + Args: + model: model to be scaled + overwrite: whether to overwrite existing scaling factors + submodel_scalers: ComponentMap of Scalers to use for sub-models + + Returns: + None + """ + if model.flow_direction is FlowDirection.forward: + x_in = model.properties.index_set().first() + elif model.flow_direction is FlowDirection.backward: + x_in = model.properties.index_set().last() + elif model.flow_direction is FlowDirection.notSet: + raise RuntimeError( + "Scaler called on ControlVolume1D without a flow " + "direction set. The unit model containing the ControlVolume1D " + "should use the add_geometry method to specify a flow direction " + "as part of model construction." + ) + else: + raise BurntToast( + "Unknown flow direction specified. This indicates " + "a new flow direction was added without support being " + "extended to the scaler. Please contact the IDAES " + "development team with this error." + ) + + self.propagate_state_scaling( + target_state=model.properties, + source_state=model.properties[x_in], + overwrite=overwrite, + ) + self.call_submodel_scaler_method( + submodel=model.properties, + submodel_scalers=submodel_scalers, + method="variable_scaling_routine", + overwrite=overwrite, + ) + for v in model.area.values(): + self.scale_variable_by_default(v, overwrite=overwrite) + + self.scale_variable_by_default(model.length, overwrite=overwrite) + if hasattr(model, "phase_fraction"): + for v in model.phase_fraction.values(): + self.scale_variable_by_default(v, overwrite=overwrite) + + super().variable_scaling_routine( + model, overwrite=overwrite, submodel_scalers=submodel_scalers + ) + + if hasattr(model, "_flow_terms"): # pylint: disable=protected-access + for idx, v in model._flow_terms.items(): + self.scale_variable_by_definition_constraint( + v, model.material_flow_linking_constraints[idx], overwrite=overwrite + ) + + if hasattr(model, "material_flow_dx"): + for idx, v in model._flow_terms.items(): # pylint: disable=protected-access + # As domain is normalized, derivative should have same + # scale as flow + self.scale_variable_by_component( + model.material_flow_dx[idx], v, overwrite=overwrite + ) + + # TODO elemental flows + # if hasattr(self, "elemental_flow_term"): + # for (t, x, e), v in self.elemental_flow_term.items(): + # flow_basis = self.properties[t, x].get_material_flow_basis() + + # sf = iscale.min_scaling_factor( + # [ + # self.properties[t, x].get_material_density_terms(p, j) + # for (p, j) in phase_component_set + # ], + # default=1, + # warning=True, + # ) + # if flow_basis == MaterialFlowBasis.molar: + # sf *= 1 + # elif flow_basis == MaterialFlowBasis.mass: + # # MW scaling factor is the inverse of its value + # sf *= value(self.properties[t, x].mw_comp[j]) + + # iscale.set_scaling_factor(v, sf) + + # if hasattr(self, "elemental_flow_dx"): + # for (t, x, e), v in self.elemental_flow_dx.items(): + # if iscale.get_scaling_factor(v) is None: + # # As domain is normalized, scale should be equal to flow + # sf = iscale.get_scaling_factor(self.elemental_flow_term[t, x, e]) + # iscale.set_scaling_factor(v, sf) + + if hasattr(model, "_enthalpy_flow"): + for ( + idx, + v, + ) in model._enthalpy_flow.items(): # pylint: disable=protected-access + self.scale_variable_by_definition_constraint( + v, model.enthalpy_flow_linking_constraint[idx], overwrite=overwrite + ) + + if hasattr(model, "enthalpy_flow_dx"): + for ( + idx, + v, + ) in model._enthalpy_flow.items(): # pylint: disable=protected-access + # Normalized domain, so scale should be the same as flow + # TODO is this correct? + self.scale_variable_by_component( + model.enthalpy_flow_dx[idx], v, overwrite=overwrite + ) + if hasattr(model, "pressure_dx"): + for (t, x), v in model.pressure_dx.items(): + self.scale_variable_by_component( + v, model.properties[t, x].pressure, overwrite=overwrite + ) + + def constraint_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None + ): + """ + Routine to apply scaling factors to constraints in model. + + Derived classes must overload this method. + + Args: + model: model to be scaled + overwrite: whether to overwrite existing scaling factors + submodel_scalers: ComponentMap of Scalers to use for sub-models + + Returns: + None + """ + self.call_submodel_scaler_method( + submodel=model.properties, + submodel_scalers=submodel_scalers, + method="constraint_scaling_routine", + overwrite=overwrite, + ) + + super().constraint_scaling_routine( + model, overwrite=overwrite, submodel_scalers=submodel_scalers + ) + + if hasattr(model, "material_flow_linking_constraints"): + for idx, v in model._flow_terms.items(): # pylint: disable=protected-access + self.scale_constraint_by_component( + model.material_flow_linking_constraints[idx], v, overwrite=overwrite + ) + + # TODO element balance + # if hasattr(model, "elemental_flow_constraint"): + # for idx, c in model.elemental_flow_constraint.items(): + # self.scale_constraint_by_component( + # c, + # model.elemental_flow_term[idx], + # overwrite=overwrite + # ) + + if hasattr(model, "enthalpy_flow_linking_constraint"): + for idx, v in model._enthalpy_flow.items(): + self.scale_constraint_by_component( + model.enthalpy_flow_linking_constraint[idx], v, overwrite=overwrite + ) + + if hasattr(model, "material_flow_dx_disc_eq"): + for idx, c in model.material_flow_dx_disc_eq.items(): + self.scale_constraint_by_component( + c, model.material_flow_dx[idx], overwrite=overwrite + ) + + if hasattr(model, "_flow_terms_length_domain_cont_eq"): + for ( + idx, + c, + ) in ( + model._flow_terms_length_domain_cont_eq.items() + ): # pylint: disable=protected-access + self.scale_constraint_by_component( + c, + model._flow_terms[idx], # pylint: disable=protected-access + overwrite=overwrite, + ) + + if hasattr(model, "enthalpy_flow_dx_disc_eq"): + for idx, c in model.enthalpy_flow_dx_disc_eq.items(): + self.scale_constraint_by_component( + c, model.enthalpy_flow_dx[idx], overwrite=overwrite + ) + + if hasattr(model, "_enthalpy_flow_length_domain_cont_eq"): + for ( + idx, + c, + ) in ( + model._enthalpy_flow_length_domain_cont_eq.items() + ): # pylint: disable=protected-access + self.scale_constraint_by_component( + c, + model._enthalpy_flow[idx], # pylint: disable=protected-access + overwrite=overwrite, + ) + + if hasattr(model, "pressure_dx_disc_eq"): + for idx, c in model.pressure_dx_disc_eq.items(): + self.scale_constraint_by_component( + c, model.pressure_dx[idx], overwrite=overwrite + ) + + if hasattr(model, "pressure_length_domain_cont_eq"): + for (t, x), c in model.pressure_length_domain_cont_eq.items(): + self.scale_constraint_by_component( + c, model.properties[t, x].pressure, overwrite=overwrite + ) + + # TODO element flow + # if hasattr(model, "element_flow_dx_disc_eq"): + # for idx, c in model.element_flow_dx_disc_eq.items(): + # self.scale_constraint_by_component( + # c, + # model.element_flow_dx[idx], + # overwrite=overwrite + # ) + # element_flow_dx_cont_eq + + @declare_process_block_class( "ControlVolume1DBlock", doc=""" @@ -99,6 +355,8 @@ class ControlVolume1DBlockData(ControlVolumeBlockData): specified in the chosen property package. """ + default_scaler = ControlVolume1DScaler + CONFIG = ControlVolumeBlockData.CONFIG() CONFIG.declare( "area_definition", @@ -155,6 +413,13 @@ class ControlVolume1DBlockData(ControlVolumeBlockData): ), ) + @property + def flow_direction(self): + """ + Property giving the flow direction in the ControlVolume1D + """ + return self._flow_direction + def build(self): """ Build method for ControlVolume1DBlock blocks. @@ -166,6 +431,7 @@ def build(self): super(ControlVolume1DBlockData, self).build() self._validate_config_args() + self._flow_direction = FlowDirection.notSet def _validate_config_args(self): # Validate DAE config arguments @@ -2437,6 +2703,8 @@ def calculate_scaling_factors(self): ) iscale.constraint_scaling_transform(c, sf, overwrite=False) + # TODO Inherent reaction stoichiometry constraint missing + if hasattr(self, "material_balances"): mb_type = self._constructed_material_balance_type if mb_type == MaterialBalanceType.componentPhase: diff --git a/idaes/core/base/control_volume_base.py b/idaes/core/base/control_volume_base.py index a1f8180550..6868364402 100644 --- a/idaes/core/base/control_volume_base.py +++ b/idaes/core/base/control_volume_base.py @@ -48,7 +48,7 @@ _log = idaeslog.getLogger(__name__) -__author__ = "Andrew Lee" +__author__ = "Andrew Lee, Douglas Allan" # Enumerate options for material balances @@ -99,6 +99,7 @@ class FlowDirection(Enum): Enum indicating direction of flow. """ + notSet = 0 forward = 1 backward = 2 @@ -110,6 +111,18 @@ class ControlVolumeScalerBase(CustomScalerBase): # TODO can we extend this to the Mixer, Separator, and MSContactor? + # String name of the reference state block to use when scaling + # balance equations and material/heat transfer terms. This should be + # set by the scaler object that inherits from this base + _state_block_ref = None + + # Attribute name to use as a weight when scaling material and energy + # terms. Presently (9/25/25), this attribute exists to take into + # account the fact that all the material and energy terms in the + # ControlVolume1D are given on the basis of material or energy per + # unit length, so we want to weigh them accordingly. + _weight_attr_name = None + def variable_scaling_routine( self, model, overwrite: bool = False, submodel_scalers=None ): @@ -126,21 +139,39 @@ def variable_scaling_routine( Returns: None """ - if hasattr(model, "properties_out"): - # ControlVolume0D - phase_list = model.properties_out.phase_list - phase_component_set = model.properties_out.phase_component_set - props = model.properties_out - elif hasattr(model, "properties"): - # ControlVolume1D - phase_list = model.properties.phase_list - phase_component_set = model.properties.phase_component_set - props = model.properties + if self._state_block_ref is None: + raise AttributeError( + "The _state_block_ref attribute was not overridden by the " + "class inheriting from ControlVolumeScalerBase." + ) else: - raise RuntimeError( - "ControlVolumeScalerBase can scale only the ControlVolume0D " - "and ControlVolume1D classes." + props = getattr(model, self._state_block_ref) + phase_list = props.phase_list + phase_component_set = props.phase_component_set + + phase_equilibrium_idx = getattr( + props.params, + "phase_equilibrium_idx", + None, # Default value if attr does not exist ) + phase_equilibrium_list = getattr( + props.params, + "phase_equilibrium_list", + None, # Default value if attr does not exist + ) + + if self._weight_attr_name is None: + weight = 1 + else: + # For ControlVolume1D, the weight is L and the + # terms have units of material per length or + # energy per length. The scaling factor of L has + # units of 1 / L, so we want to divide by the scaling + # factor to render the material or energy terms + # dimensionless. + weight_attr = getattr(model, self._weight_attr_name) + weight = 1 / self.get_scaling_factor(weight_attr, default=1, warning=True) + idx0 = props.index_set().first() params = props[idx0].params if hasattr(model, "reactions"): @@ -160,6 +191,7 @@ def variable_scaling_routine( ) # Material accumulation should be scaled by a global method for scaling + # time DerivativeVars if hasattr(model, "material_accumulation"): pass @@ -183,7 +215,7 @@ def variable_scaling_routine( props[prop_idx].get_material_flow_terms(p, j) ) self.set_component_scaling_factor( - rate_rxn_gen[idx], 1 / nom, overwrite=overwrite + rate_rxn_gen[idx], weight / nom, overwrite=overwrite ) # Extent of reaction scaling is based on the species in the @@ -208,7 +240,9 @@ def variable_scaling_routine( f"Reaction {rxn} has no nonzero stoichiometric coefficient." ) # Note this scaling works only if we don't - # have multiple reactions cancelling each other out + # have multiple reactions cancelling each other out. + # No need to weight here because the rate_rxn_gen + # scaling factor is already weighted self.set_component_scaling_factor( model.rate_reaction_extent[prop_idx, rxn], 1 / nom_rxn, @@ -230,7 +264,7 @@ def variable_scaling_routine( props[prop_idx].get_material_flow_terms(p, j) ) self.set_component_scaling_factor( - equil_rxn_gen[idx], 1 / nom, overwrite=overwrite + equil_rxn_gen[idx], weight / nom, overwrite=overwrite ) # Extent of reaction scaling is based on the species in the @@ -256,6 +290,8 @@ def variable_scaling_routine( ) # Note this scaling works only if we don't # have multiple reactions cancelling each other out + # No need to weight here because the equil_rxn_gen + # scaling factor is already weighted self.set_component_scaling_factor( model.equilibrium_reaction_extent[prop_idx, rxn], 1 / nom_rxn, @@ -277,7 +313,7 @@ def variable_scaling_routine( props[prop_idx].get_material_flow_terms(p, j) ) self.set_component_scaling_factor( - inh_rxn_gen[idx], 1 / nom, overwrite=overwrite + inh_rxn_gen[idx], weight / nom, overwrite=overwrite ) # Extent of reaction scaling is based on the species in the @@ -303,6 +339,8 @@ def variable_scaling_routine( ) # Note this scaling works only if we don't # have multiple reactions cancelling each other out + # No need to weight here because the inh_rxn_gen + # scaling factor is already weighted self.set_component_scaling_factor( model.inherent_reaction_extent[prop_idx, rxn], 1 / nom_rxn, @@ -317,10 +355,27 @@ def variable_scaling_routine( ) self.set_component_scaling_factor( model.mass_transfer_term[prop_idx, p, j], - 1 / nom, + weight / nom, overwrite=overwrite, ) + if hasattr(model, "phase_equilibrium_generation"): + for prop_idx in props: + for pe_idx in phase_equilibrium_idx: + j, pp = phase_equilibrium_list[pe_idx] + nom1 = self.get_expression_nominal_value( + props[prop_idx].get_material_flow_terms(pp[0], j) + ) + nom2 = self.get_expression_nominal_value( + props[prop_idx].get_material_flow_terms(pp[1], j) + ) + nom = min(nom1, nom2) + self.set_component_scaling_factor( + model.phase_equilibrium_generation[prop_idx, pe_idx], + 1 / nom, + overwrite=False, + ) + # Set scaling factors for element balance variables # TODO # if hasattr(model, "elemental_flow_out"): @@ -415,25 +470,27 @@ def variable_scaling_routine( nom = max(nom_list) if hasattr(model, "heat"): self.set_component_scaling_factor( - model.heat[prop_idx], 1 / nom, overwrite=overwrite + model.heat[prop_idx], weight / nom, overwrite=overwrite ) if hasattr(model, "work"): self.set_component_scaling_factor( - model.work[prop_idx], 1 / nom, overwrite=overwrite + model.work[prop_idx], weight / nom, overwrite=overwrite ) if hasattr(model, "enthalpy_transfer"): self.set_component_scaling_factor( - model.enthalpy_transfer[prop_idx], 1 / nom, overwrite=overwrite + model.enthalpy_transfer[prop_idx], + weight / nom, + overwrite=overwrite, ) # Set scaling for momentum balance variables if hasattr(model, "deltaP"): for prop_idx in props: - sf_P = get_scaling_factor(props[prop_idx].pressure) - # TODO raise error if pressure scaling factor - # isn't set + sf_P = get_scaling_factor( + props[prop_idx].pressure, default=1e-5, warning=True + ) self.set_component_scaling_factor( - model.deltaP[prop_idx], sf_P, overwrite=overwrite + model.deltaP[prop_idx], weight * sf_P, overwrite=overwrite ) def constraint_scaling_routine( @@ -452,19 +509,15 @@ def constraint_scaling_routine( Returns: None """ - if hasattr(model, "properties_out"): - # ControlVolume0D - phase_list = model.properties_out.phase_list - props = model.properties_out - elif hasattr(model, "properties"): - # ControlVolume1D - phase_list = model.properties.phase_list - props = model.properties - else: - raise RuntimeError( - "ControlVolumeScalerBase can scale only the ControlVolume0D " - "and ControlVolume1D classes." + if self._state_block_ref is None: + raise AttributeError( + "The _state_block_ref attribute was not overridden by the " + "class inheriting from ControlVolumeScalerBase." ) + else: + props = getattr(model, self._state_block_ref) + phase_list = props.phase_list + pc_set = props.phase_component_set if hasattr(model, "reactions"): self.call_submodel_scaler_method( @@ -473,7 +526,7 @@ def constraint_scaling_routine( method="constraint_scaling_routine", overwrite=overwrite, ) - # Transform constraints in order of appearance + if hasattr(model, "material_holdup_calculation"): for idx in model.material_holdup_calculation: self.scale_constraint_by_component( @@ -498,30 +551,76 @@ def constraint_scaling_routine( overwrite=overwrite, ) + inh_rxn_con = None if hasattr(model, "inherent_reaction_stoichiometry_constraint"): - for idx in model.inherent_reaction_stoichiometry_constraint: + # ControlVolume0D and ControlVolume1D + inh_rxn_con = model.inherent_reaction_stoichiometry_constraint + elif hasattr(model, "inherent_reaction_constraint"): + # Mixer + inh_rxn_con = model.inherent_reaction_constraint + + if inh_rxn_con is not None: + for idx in inh_rxn_con: self.scale_constraint_by_component( - model.inherent_reaction_stoichiometry_constraint[idx], + inh_rxn_con[idx], model.inherent_reaction_generation[idx], overwrite=overwrite, ) + mb_eqn = None if hasattr(model, "material_balances"): + # ControlVolume0D and ControlVolume1D + mb_eqn = model.material_balances + elif hasattr(model, "material_mixing_equations"): + # Mixer + mb_eqn = model.material_mixing_equations + + if mb_eqn is not None: mb_type = model._constructed_material_balance_type # pylint: disable=W0212 - if ( - mb_type == MaterialBalanceType.componentPhase - or mb_type == MaterialBalanceType.componentTotal - ): - for idx in model.material_balances: - self.scale_constraint_by_nominal_value( - model.material_balances[idx], - scheme=ConstraintScalingScheme.inverseMaximum, - overwrite=overwrite, + if mb_type == MaterialBalanceType.componentTotal: + for idx in mb_eqn: + c = idx[-1] + nom_list = [] + for p in phase_list: + nom_list.append( + self.get_expression_nominal_value( + props[idx[:-1]].get_material_flow_terms(p, c) + ) + ) + nom = max(nom_list) + self.set_component_scaling_factor( + mb_eqn[idx], 1 / nom, overwrite=overwrite + ) + elif mb_type == MaterialBalanceType.componentPhase: + for idx in mb_eqn: + p = idx[-2] + c = idx[-1] + nom = self.get_expression_nominal_value( + props[idx[:-2]].get_material_flow_terms(p, c) + ) + self.set_component_scaling_factor( + mb_eqn[idx], 1 / nom, overwrite=overwrite + ) + elif mb_type == MaterialBalanceType.total: + for idx in mb_eqn: + nom_list = [] + for p, c in pc_set: + nom_list.append( + self.get_expression_nominal_value( + props[idx[:-1]].get_material_flow_terms(p, c) + ) + ) + nom = max(nom_list) + self.set_component_scaling_factor( + mb_eqn[idx], 1 / nom, overwrite=overwrite ) else: # There are some other material balance types but they create # constraints with different names. - _log.warning(f"Unknown material balance type {mb_type}") + _log.warning( + f"Unknown material balance type {mb_type}. It cannot be " + "automatically scaled." + ) # TODO element balances # if hasattr(self, "element_balances"): @@ -536,8 +635,16 @@ def constraint_scaling_routine( # sf = iscale.get_scaling_factor(self.element_holdup[t, e]) # iscale.constraint_scaling_transform(c, sf, overwrite=False) + eb_eqn = None if hasattr(model, "enthalpy_balances"): - for idx in props: + eb_eqn = model.enthalpy_balances + elif hasattr(model, "enthalpy_mixing_equations"): + eb_eqn = model.enthalpy_mixing_equations + + if eb_eqn is not None: + # Phase enthalpy balances are not implemented + # as of 9/26/25 + for idx in eb_eqn: nom_list = [] for p in phase_list: nom_list.append( @@ -547,7 +654,7 @@ def constraint_scaling_routine( ) nom = max(nom_list) self.set_component_scaling_factor( - model.enthalpy_balances[idx], 1 / nom, overwrite=overwrite + eb_eqn[idx], 1 / nom, overwrite=overwrite ) if hasattr(model, "energy_holdup_calculation"): @@ -558,11 +665,16 @@ def constraint_scaling_routine( overwrite=overwrite, ) + pb_eqn = None if hasattr(model, "pressure_balance"): - for con in model.pressure_balance.values(): - self.scale_constraint_by_nominal_value( + # ControlVolume0D and ControlVolume1D + pb_eqn = model.pressure_balance + + if pb_eqn is not None: + for idx, con in pb_eqn.items(): + self.scale_constraint_by_component( con, - scheme=ConstraintScalingScheme.inverseMaximum, + props[idx].pressure, overwrite=overwrite, ) diff --git a/idaes/core/base/tests/test_control_volume_0d_scaling_legacy.py b/idaes/core/base/tests/test_control_volume_0d_legacy_scaling.py similarity index 100% rename from idaes/core/base/tests/test_control_volume_0d_scaling_legacy.py rename to idaes/core/base/tests/test_control_volume_0d_legacy_scaling.py diff --git a/idaes/core/base/tests/test_control_volume_0d_scaler_object.py b/idaes/core/base/tests/test_control_volume_0d_scaler_object.py index a0f7831b34..11afc0e9ba 100644 --- a/idaes/core/base/tests/test_control_volume_0d_scaler_object.py +++ b/idaes/core/base/tests/test_control_volume_0d_scaler_object.py @@ -34,6 +34,10 @@ from idaes.core.base.control_volume0d import ControlVolume0DScaler +def approx(x): + return pytest.approx(x, rel=1e-15, abs=0) + + @pytest.mark.unit def test_basic_scaling(): m = pyo.ConcreteModel() @@ -179,7 +183,7 @@ def test_full_auto_scaling_dynamic(): assert get_scaling_factor(v) == 1 / 71 for v in cv.material_holdup.values(): - assert get_scaling_factor(v) == pytest.approx(10 / (47 * 71), rel=1e-15, abs=0) + assert get_scaling_factor(v) == approx(10 / (47 * 71)) for v in cv.rate_reaction_generation.values(): assert get_scaling_factor(v) == 1 / 43 @@ -200,7 +204,7 @@ def test_full_auto_scaling_dynamic(): assert get_scaling_factor(v) == 10 for v in cv.energy_holdup.values(): - assert get_scaling_factor(v) == pytest.approx(10 / (41 * 71), rel=1e-15, abs=0) + assert get_scaling_factor(v) == approx(10 / (41 * 71)) for v in cv.heat.values(): assert get_scaling_factor(v) == 1 / 37 @@ -218,7 +222,7 @@ def test_full_auto_scaling_dynamic(): assert get_scaling_factor(c) == 1 for c in cv.material_holdup_calculation.values(): - assert get_scaling_factor(c) == pytest.approx(10 / (47 * 71), rel=1e-15, abs=0) + assert get_scaling_factor(c) == approx(10 / (47 * 71)) for c in cv.rate_reaction_stoichiometry_constraint.values(): assert get_scaling_factor(c) == 1 / 43 @@ -233,7 +237,7 @@ def test_full_auto_scaling_dynamic(): assert get_scaling_factor(c) == 1 / 37 for c in cv.energy_holdup_calculation.values(): - assert get_scaling_factor(c) == pytest.approx(10 / (41 * 71), rel=1e-15, abs=0) + assert get_scaling_factor(c) == approx(10 / (41 * 71)) for c in cv.pressure_balance.values(): assert get_scaling_factor(c) == 1 / 13 @@ -245,6 +249,101 @@ def test_full_auto_scaling_dynamic(): assert len(list_unscaled_constraints(m)) == 18 +@pytest.mark.unit +def test_full_auto_scaling_mbtype_phase(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp) + m.fs.cv = ControlVolume0DBlock( + property_package=m.fs.pp, + reaction_package=m.fs.rp, + ) + m.fs.cv.add_geometry() + m.fs.cv.add_state_blocks(has_phase_equilibrium=True) + m.fs.cv.add_reaction_blocks(has_equilibrium=True) + + m.fs.cv.add_material_balances( + balance_type=MaterialBalanceType.componentPhase, + has_rate_reactions=True, + has_equilibrium_reactions=True, + has_phase_equilibrium=True, + has_mass_transfer=True, + ) + + m.fs.cv.add_energy_balances( + balance_type=EnergyBalanceType.enthalpyTotal, + has_heat_of_reaction=True, + has_heat_transfer=True, + has_work_transfer=True, + has_enthalpy_transfer=True, + ) + + m.fs.cv.add_momentum_balances( + balance_type=MomentumBalanceType.pressureTotal, has_pressure_change=True + ) + + scaler_obj = ControlVolume0DScaler() + scaler_obj.default_scaling_factors["volume"] = 1 / 71 + + scaler_obj.scale_model(m.fs.cv) + + assert len(list_unscaled_variables(m, include_fixed=True)) == 0 + assert len(list_unscaled_constraints(m)) == 0 + + cv = m.fs.cv + + # Variables + for v in cv.volume.values(): + assert get_scaling_factor(v) == 1 / 71 + + for v in cv.rate_reaction_generation.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.rate_reaction_extent.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.equilibrium_reaction_generation.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.equilibrium_reaction_extent.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.mass_transfer_term.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.heat.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.work.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.enthalpy_transfer.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.deltaP.values(): + assert get_scaling_factor(v) == 1 / 13 + + # Constraints + for c in cv.rate_reaction_stoichiometry_constraint.values(): + assert get_scaling_factor(c) == 1 / 43 + + for c in cv.equilibrium_reaction_stoichiometry_constraint.values(): + assert get_scaling_factor(c) == 1 / 43 + + for c in cv.material_balances.values(): + assert get_scaling_factor(c) == 1 / 43 + + for c in cv.enthalpy_balances.values(): + assert get_scaling_factor(c) == 1 / 37 + + for c in cv.pressure_balance.values(): + assert get_scaling_factor(c) == 1 / 13 + + for c in cv.phase_equilibrium_generation.values(): + assert get_scaling_factor(c) == 1 / 43 + + # TODO test element balance when implemented # @pytest.mark.unit # def test_full_auto_scaling_mbtype_element(): diff --git a/idaes/core/base/tests/test_control_volume_1d.py b/idaes/core/base/tests/test_control_volume_1d.py index 073ab5cce0..4475dd826e 100644 --- a/idaes/core/base/tests/test_control_volume_1d.py +++ b/idaes/core/base/tests/test_control_volume_1d.py @@ -195,6 +195,10 @@ def test_add_geometry_default(): finite_elements=10, ) + assert m.fs.cv.flow_direction == FlowDirection.notSet + with pytest.raises(AttributeError): + m.fs.cv.flow_direction = FlowDirection.backward + m.fs.cv.add_geometry() assert isinstance(m.fs.cv.length_domain, ContinuousSet) @@ -206,6 +210,9 @@ def test_add_geometry_default(): assert len(m.fs.cv.length) == 1.0 assert m.fs.cv.length.value == 1.0 assert m.fs.cv._flow_direction == FlowDirection.forward + assert m.fs.cv.flow_direction == FlowDirection.forward + with pytest.raises(AttributeError): + m.fs.cv.flow_direction = FlowDirection.backward @pytest.mark.unit diff --git a/idaes/core/base/tests/test_control_volume_1d_scaling.py b/idaes/core/base/tests/test_control_volume_1d_legacy_scaling.py similarity index 100% rename from idaes/core/base/tests/test_control_volume_1d_scaling.py rename to idaes/core/base/tests/test_control_volume_1d_legacy_scaling.py diff --git a/idaes/core/base/tests/test_control_volume_1d_scaler_object.py b/idaes/core/base/tests/test_control_volume_1d_scaler_object.py new file mode 100644 index 0000000000..685136bc76 --- /dev/null +++ b/idaes/core/base/tests/test_control_volume_1d_scaler_object.py @@ -0,0 +1,909 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2024 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Tests for ControlVolume1D scaler object. + +Author: Douglas Allan +""" +import pytest +import re +import pyomo.environ as pyo +from idaes.core import ( + ControlVolume1DBlock, + MaterialBalanceType, + EnergyBalanceType, + MomentumBalanceType, + FlowsheetBlock, + FlowDirection, +) +from idaes.core.base.control_volume1d import ControlVolume1DScaler +from idaes.core.util.testing import ( + PhysicalParameterTestBlock, + ReactionParameterTestBlock, +) +from idaes.core.util.exceptions import BurntToast +from idaes.core.scaling import get_scaling_factor +from idaes.core.scaling.util import list_unscaled_variables, list_unscaled_constraints + + +def approx(x): + return pytest.approx(x, rel=1e-15, abs=0) + + +# ----------------------------------------------------------------------------- +# Basic tests +@pytest.mark.unit +def test_basic_scaling(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + # Set flag to include inherent reactions + m.fs.pp._has_inherent_reactions = True + + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + + m.fs.cv.add_geometry() + m.fs.cv.add_state_blocks(has_phase_equilibrium=False) + + m.fs.cv.add_material_balances( + balance_type=MaterialBalanceType.componentTotal, has_phase_equilibrium=False + ) + + m.fs.cv.add_energy_balances(balance_type=EnergyBalanceType.enthalpyTotal) + + m.fs.cv.add_momentum_balances( + balance_type=MomentumBalanceType.pressureTotal, has_pressure_change=True + ) + + m.fs.cv.apply_transformation() + + assert m.fs.cv.default_scaler is ControlVolume1DScaler + + scaler_obj = ControlVolume1DScaler() + scaler_obj.default_scaling_factors["area"] = 1 / 83 + scaler_obj.default_scaling_factors["length"] = 1 / 599 + scaler_obj.scale_model(m.fs.cv) + + # check scaling on select variables + assert get_scaling_factor(m.fs.cv.area) == 1 / 83 + + for v in m.fs.cv._flow_terms.values(): + assert get_scaling_factor(v) == approx(1 / 43) + + for v in m.fs.cv.material_flow_dx.values(): + assert get_scaling_factor(v) == approx(1 / 43) + + for v in m.fs.cv._enthalpy_flow.values(): + assert get_scaling_factor(v) == approx(1 / 37) + + for v in m.fs.cv.enthalpy_flow_dx.values(): + assert get_scaling_factor(v) == approx(1 / 37) + + for v in m.fs.cv.deltaP.values(): + assert get_scaling_factor(v) == approx(599 / 13) + for v in m.fs.cv.pressure_dx.values(): + get_scaling_factor(v) == approx(1 / 13) + + for t in m.fs.time: + for x in m.fs.cv.length_domain: + assert get_scaling_factor(m.fs.cv.properties[t, x].flow_vol) == 1 / 3 + + # check scaling on mass, energy, and pressure balances. + for c in m.fs.cv.material_balances.values(): + assert get_scaling_factor(c) == 1 / 43 + for c in m.fs.cv.enthalpy_balances.values(): + # this uses the minimum enthalpy flow term scale + assert get_scaling_factor(c) == 1 / 37 + for c in m.fs.cv.pressure_balance.values(): + # This uses the inlet pressure scale + assert get_scaling_factor(c) == 1 / 13 + + +@pytest.mark.unit +def test_user_set_scaling_FD_forward(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + m.fs.cv.add_geometry() + m.fs.cv.add_state_blocks(has_phase_equilibrium=False) + m.fs.cv.add_material_balances( + balance_type=MaterialBalanceType.componentTotal, has_phase_equilibrium=False + ) + m.fs.cv.add_energy_balances( + balance_type=EnergyBalanceType.enthalpyTotal, + has_heat_transfer=True, + has_work_transfer=True, + ) + # add momentum balance + m.fs.cv.add_momentum_balances( + balance_type=MomentumBalanceType.pressureTotal, has_pressure_change=True + ) + + m.fs.cv.apply_transformation() + + scaler_obj = ControlVolume1DScaler() + + # Scaling factors should get propagated from inlet state block to the other + # state blocks + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].temperature, 1 / 1009 + ) + scaler_obj.set_component_scaling_factor(m.fs.cv.properties[0, 0].pressure, 1 / 1301) + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].flow_mol_phase_comp["p1", "c1"], 1 / 1117 + ) + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].flow_mol_phase_comp["p2", "c1"], 1 / 1931 + ) + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].flow_mol_phase_comp["p1", "c2"], 1 / 1549 + ) + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].flow_mol_phase_comp["p2", "c2"], 1 / 1999 + ) + + scaler_obj.set_component_scaling_factor(m.fs.cv.properties[0, 1].pressure, 1 / 2251) + + scaler_obj.default_scaling_factors["area"] = 1 / 83 + scaler_obj.default_scaling_factors["length"] = 1 / 599 + + # The scaling factors used for this test were selected to be easy values to + # test, they do not represent typical scaling factors. + for v in m.fs.cv.heat.values(): + scaler_obj.set_component_scaling_factor(v, 1 / 73) + for v in m.fs.cv.work.values(): + scaler_obj.set_component_scaling_factor(v, 1 / 79) + scaler_obj.set_component_scaling_factor(m.fs.cv.heat[0, 0], 17, overwrite=True) + + scaler_obj.scale_model(m.fs.cv) + for t in m.fs.time: + for x in m.fs.cv.length_domain: + if t == 0 and x == 0: + assert get_scaling_factor(m.fs.cv.heat[t, x]) == 17 + else: + assert get_scaling_factor(m.fs.cv.heat[t, x]) == 1 / 73 + assert get_scaling_factor(m.fs.cv.work[t, x]) == 1 / 79 + + assert get_scaling_factor(m.fs.cv.properties[t, x].temperature) == 1 / 1009 + if t == 0 and x == 1: + # Propagation from inlet state block shouldn't overwrite existing scaling factor + assert get_scaling_factor(m.fs.cv.properties[t, x].pressure) == 1 / 2251 + else: + assert get_scaling_factor(m.fs.cv.properties[t, x].pressure) == 1 / 1301 + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p1", "c1"] + ) + == 1 / 1117 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p2", "c1"] + ) + == 1 / 1931 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p1", "c2"] + ) + == 1 / 1549 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p2", "c2"] + ) + == 1 / 1999 + ) + + +@pytest.mark.unit +def test_user_set_scaling_FD_backward(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + m.fs.cv.add_geometry(flow_direction=FlowDirection.backward) + m.fs.cv.add_state_blocks(has_phase_equilibrium=False) + m.fs.cv.add_material_balances( + balance_type=MaterialBalanceType.componentTotal, has_phase_equilibrium=False + ) + m.fs.cv.add_energy_balances( + balance_type=EnergyBalanceType.enthalpyTotal, + has_heat_transfer=True, + has_work_transfer=True, + ) + # add momentum balance + m.fs.cv.add_momentum_balances( + balance_type=MomentumBalanceType.pressureTotal, has_pressure_change=True + ) + + m.fs.cv.apply_transformation() + + scaler_obj = ControlVolume1DScaler() + + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].temperature, 1 / 1009 + ) + scaler_obj.set_component_scaling_factor(m.fs.cv.properties[0, 0].pressure, 1 / 1301) + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].flow_mol_phase_comp["p1", "c1"], 1 / 1117 + ) + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].flow_mol_phase_comp["p2", "c1"], 1 / 1931 + ) + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].flow_mol_phase_comp["p1", "c2"], 1 / 1549 + ) + scaler_obj.set_component_scaling_factor( + m.fs.cv.properties[0, 0].flow_mol_phase_comp["p2", "c2"], 1 / 1999 + ) + + scaler_obj.set_component_scaling_factor(m.fs.cv.properties[0, 1].pressure, 1 / 2251) + + scaler_obj.default_scaling_factors["area"] = 1 / 83 + scaler_obj.default_scaling_factors["length"] = 1 / 599 + + # The scaling factors used for this test were selected to be easy values to + # test, they do not represent typical scaling factors. + for v in m.fs.cv.heat.values(): + scaler_obj.set_component_scaling_factor(v, 1 / 73) + for v in m.fs.cv.work.values(): + scaler_obj.set_component_scaling_factor(v, 1 / 79) + scaler_obj.set_component_scaling_factor(m.fs.cv.heat[0, 0], 17, overwrite=True) + + scaler_obj.scale_model(m.fs.cv) + # Since the inlet state block is [0, 1] and the only scaling factor set on + # that block is for pressure, only that scaling factor is propagated to the + # other state blocks. All other state variables use default values + for t in m.fs.time: + for x in m.fs.cv.length_domain: + if t == 0 and x == 0: + assert get_scaling_factor(m.fs.cv.heat[t, x]) == 17 + else: + assert get_scaling_factor(m.fs.cv.heat[t, x]) == 1 / 73 + assert get_scaling_factor(m.fs.cv.work[t, x]) == 1 / 79 + + if t == 0 and x == 0: + # Scaling factors were specifically set at [0, 0], so those + # shouldn't get overwritten + assert ( + get_scaling_factor(m.fs.cv.properties[t, x].temperature) == 1 / 1009 + ) + assert get_scaling_factor(m.fs.cv.properties[t, x].pressure) == 1 / 1301 + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p1", "c1"] + ) + == 1 / 1117 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p2", "c1"] + ) + == 1 / 1931 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p1", "c2"] + ) + == 1 / 1549 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p2", "c2"] + ) + == 1 / 1999 + ) + else: + assert ( + get_scaling_factor(m.fs.cv.properties[t, x].temperature) == 1 / 17 + ) + assert get_scaling_factor(m.fs.cv.properties[t, x].pressure) == 1 / 2251 + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p1", "c1"] + ) + == 1 / 7 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p2", "c1"] + ) + == 1 / 7 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p1", "c2"] + ) + == 1 / 7 + ) + assert ( + get_scaling_factor( + m.fs.cv.properties[t, x].flow_mol_phase_comp["p2", "c2"] + ) + == 1 / 7 + ) + + +@pytest.mark.unit +def test_user_set_scaling_FD_not_set(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + + scaler_obj = ControlVolume1DScaler() + + with pytest.raises( + RuntimeError, + match=re.escape( + "Scaler called on ControlVolume1D without a flow " + "direction set. The unit model containing the ControlVolume1D " + "should use the add_geometry method to specify a flow direction " + "as part of model construction." + ), + ): + scaler_obj.scale_model(m.fs.cv) + + +@pytest.mark.unit +def test_user_set_scaling_FD_burnt_toast(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + m.fs.cv._flow_direction = "foo" + + scaler_obj = ControlVolume1DScaler() + + with pytest.raises( + BurntToast, + match=re.escape( + "Unknown flow direction specified. This indicates " + "a new flow direction was added without support being " + "extended to the scaler. Please contact the IDAES " + "development team with this error." + ), + ): + scaler_obj.scale_model(m.fs.cv) + + +@pytest.mark.unit +def test_full_auto_scaling(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp) + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + reaction_package=m.fs.rp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + m.fs.cv.add_geometry() + m.fs.cv.add_state_blocks(has_phase_equilibrium=True) + m.fs.cv.add_reaction_blocks(has_equilibrium=True) + + m.fs.cv.add_material_balances( + balance_type=MaterialBalanceType.componentTotal, + has_rate_reactions=True, + has_equilibrium_reactions=True, + has_phase_equilibrium=True, + has_mass_transfer=True, + ) + + m.fs.cv.add_energy_balances( + balance_type=EnergyBalanceType.enthalpyTotal, + has_heat_of_reaction=True, + has_heat_transfer=True, + has_work_transfer=True, + has_enthalpy_transfer=True, + ) + + m.fs.cv.add_momentum_balances( + balance_type=MomentumBalanceType.pressureTotal, has_pressure_change=True + ) + + m.fs.cv.apply_transformation() + + scaler_obj = ControlVolume1DScaler() + scaler_obj.default_scaling_factors["area"] = 1 / 83 + scaler_obj.default_scaling_factors["length"] = 1 / 599 + scaler_obj.scale_model(m.fs.cv) + + cv = m.fs.cv + + # Variables + assert get_scaling_factor(cv.length) == 1 / 599 + for v in cv.area.values(): + assert get_scaling_factor(v) == 1 / 83 + + for v in cv._flow_terms.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv._enthalpy_flow.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.material_flow_dx.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.enthalpy_flow_dx.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.pressure_dx.values(): + assert get_scaling_factor(v) == 1 / 13 + + for v in cv.rate_reaction_generation.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.rate_reaction_extent.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.equilibrium_reaction_generation.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.equilibrium_reaction_extent.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.mass_transfer_term.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.heat.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.work.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.enthalpy_transfer.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.deltaP.values(): + assert get_scaling_factor(v) == approx(599 / 13) + + # Constraints + for v in cv.material_flow_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.enthalpy_flow_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.pressure_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 13 + + for c in cv.rate_reaction_stoichiometry_constraint.values(): + assert get_scaling_factor(c) == approx(599 / 43) + + for c in cv.equilibrium_reaction_stoichiometry_constraint.values(): + assert get_scaling_factor(c) == approx(599 / 43) + + for c in cv.material_balances.values(): + assert get_scaling_factor(c) == approx(1 / 43) + + for c in cv.enthalpy_balances.values(): + assert get_scaling_factor(c) == approx(1 / 37) + + for c in cv.pressure_balance.values(): + assert get_scaling_factor(c) == approx(1 / 13) + + assert len(list_unscaled_variables(m, include_fixed=True)) == 0 + assert len(list_unscaled_constraints(m)) == 0 + + +@pytest.mark.unit +def test_full_auto_scaling_dynamic(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=True, time_set=[0, 1], time_units=pyo.units.s) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp) + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + reaction_package=m.fs.rp, + transformation_method="dae.collocation", + transformation_scheme="LAGRANGE-LEGENDRE", + finite_elements=4, + collocation_points=3, + ) + m.fs.cv.add_geometry() + m.fs.cv.add_state_blocks(has_phase_equilibrium=True) + m.fs.cv.add_reaction_blocks(has_equilibrium=True) + + m.fs.cv.add_material_balances( + balance_type=MaterialBalanceType.componentTotal, + has_rate_reactions=True, + has_equilibrium_reactions=True, + has_phase_equilibrium=True, + has_mass_transfer=True, + ) + + m.fs.cv.add_energy_balances( + balance_type=EnergyBalanceType.enthalpyTotal, + has_heat_of_reaction=True, + has_heat_transfer=True, + has_work_transfer=True, + has_enthalpy_transfer=True, + ) + + m.fs.cv.add_momentum_balances( + balance_type=MomentumBalanceType.pressureTotal, has_pressure_change=True + ) + + m.fs.cv.apply_transformation() + m.discretizer = pyo.TransformationFactory("dae.finite_difference") + m.discretizer.apply_to(m, nfe=3, wrt=m.fs.time, scheme="BACKWARD") + + scaler_obj = ControlVolume1DScaler() + scaler_obj.default_scaling_factors["area"] = 1 / 83 + scaler_obj.default_scaling_factors["length"] = 1 / 599 + + scaler_obj.scale_model(m.fs.cv) + + # Four finite elements with three collocation points on the element's + # interior plus two boundary points (which overlap with other finite elements) + # That results in 4 points per finite element plus one point at the end of the + # final finite element. + n_points_colloc = 4 * (3 + 1) + 1 + # Unscaled variables are material accumulation and energy accumulation DerivativeVars + # Two phases with two components. Energy accumulation is also indexed by phase + + assert len(list_unscaled_variables(m, include_fixed=True)) == 4 * ( + 2 * 2 * n_points_colloc + 2 * n_points_colloc + ) + # There are variables for all 4 timepoints but only constraints for 3 timepoints + assert len(list_unscaled_constraints(m)) == 3 * ( + 2 * 2 * n_points_colloc + 2 * n_points_colloc + ) + + +@pytest.mark.unit +def test_full_auto_scaling_mbtype_phase(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp) + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + reaction_package=m.fs.rp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + m.fs.cv.add_geometry() + m.fs.cv.add_state_blocks(has_phase_equilibrium=True) + m.fs.cv.add_reaction_blocks(has_equilibrium=True) + + m.fs.cv.add_material_balances( + balance_type=MaterialBalanceType.componentPhase, + has_rate_reactions=True, + has_equilibrium_reactions=True, + has_phase_equilibrium=True, + has_mass_transfer=True, + ) + + m.fs.cv.add_energy_balances( + balance_type=EnergyBalanceType.enthalpyTotal, + has_heat_of_reaction=True, + has_heat_transfer=True, + has_work_transfer=True, + has_enthalpy_transfer=True, + ) + + m.fs.cv.add_momentum_balances( + balance_type=MomentumBalanceType.pressureTotal, has_pressure_change=True + ) + + m.fs.cv.apply_transformation() + + scaler_obj = ControlVolume1DScaler() + scaler_obj.default_scaling_factors["area"] = 1 / 83 + scaler_obj.default_scaling_factors["length"] = 1 / 599 + + scaler_obj.scale_model(m.fs.cv) + + assert len(list_unscaled_variables(m, include_fixed=True)) == 0 + assert len(list_unscaled_constraints(m)) == 0 + + cv = m.fs.cv + + # Variables + assert get_scaling_factor(cv.length) == 1 / 599 + for v in cv.area.values(): + assert get_scaling_factor(v) == 1 / 83 + + for v in cv._flow_terms.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv._enthalpy_flow.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.material_flow_dx.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.enthalpy_flow_dx.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.pressure_dx.values(): + assert get_scaling_factor(v) == 1 / 13 + + for v in cv.rate_reaction_generation.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.rate_reaction_extent.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.equilibrium_reaction_generation.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.equilibrium_reaction_extent.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.mass_transfer_term.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.heat.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.work.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.enthalpy_transfer.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.deltaP.values(): + assert get_scaling_factor(v) == approx(599 / 13) + + # Constraints + for v in cv.material_flow_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.enthalpy_flow_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.pressure_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 13 + + for c in cv.rate_reaction_stoichiometry_constraint.values(): + assert get_scaling_factor(c) == approx(599 / 43) + + for c in cv.equilibrium_reaction_stoichiometry_constraint.values(): + assert get_scaling_factor(c) == approx(599 / 43) + + for c in cv.material_balances.values(): + assert get_scaling_factor(c) == approx(1 / 43) + + for c in cv.enthalpy_balances.values(): + assert get_scaling_factor(c) == approx(1 / 37) + + for c in cv.pressure_balance.values(): + assert get_scaling_factor(c) == approx(1 / 13) + + for c in cv.phase_equilibrium_generation.values(): + assert get_scaling_factor(c) == approx(1 / 43) + + +@pytest.mark.unit +def test_auto_scaling_just_material_balance(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp) + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + reaction_package=m.fs.rp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + m.fs.cv.add_geometry() + m.fs.cv.add_state_blocks(has_phase_equilibrium=True) + m.fs.cv.add_reaction_blocks(has_equilibrium=True) + + m.fs.cv.add_material_balances( + balance_type=MaterialBalanceType.componentTotal, + has_rate_reactions=True, + has_equilibrium_reactions=True, + has_phase_equilibrium=True, + has_mass_transfer=True, + ) + + m.fs.cv.apply_transformation() + + scaler_obj = ControlVolume1DScaler() + scaler_obj.default_scaling_factors["area"] = 1 / 83 + scaler_obj.default_scaling_factors["length"] = 1 / 599 + scaler_obj.scale_model(m.fs.cv) + + cv = m.fs.cv + + # Variables + assert get_scaling_factor(cv.length) == 1 / 599 + for v in cv.area.values(): + assert get_scaling_factor(v) == 1 / 83 + + for v in cv._flow_terms.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.material_flow_dx.values(): + assert get_scaling_factor(v) == 1 / 43 + + for v in cv.rate_reaction_generation.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.rate_reaction_extent.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.equilibrium_reaction_generation.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.equilibrium_reaction_extent.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + for v in cv.mass_transfer_term.values(): + assert get_scaling_factor(v) == approx(599 / 43) + + # Constraints + for v in cv.material_flow_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 43 + for c in cv.rate_reaction_stoichiometry_constraint.values(): + assert get_scaling_factor(c) == approx(599 / 43) + + for c in cv.equilibrium_reaction_stoichiometry_constraint.values(): + assert get_scaling_factor(c) == approx(599 / 43) + + for c in cv.material_balances.values(): + assert get_scaling_factor(c) == approx(1 / 43) + + assert len(list_unscaled_variables(m, include_fixed=True)) == 0 + assert len(list_unscaled_constraints(m)) == 0 + + +@pytest.mark.unit +def test_auto_scaling_enthalpy_and_pressure_balance(): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp) + m.fs.cv = ControlVolume1DBlock( + property_package=m.fs.pp, + reaction_package=m.fs.rp, + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + finite_elements=10, + ) + m.fs.cv.add_geometry() + m.fs.cv.add_state_blocks(has_phase_equilibrium=True) + + m.fs.cv.add_energy_balances( + balance_type=EnergyBalanceType.enthalpyTotal, + has_heat_of_reaction=False, + has_heat_transfer=True, + has_work_transfer=True, + has_enthalpy_transfer=True, + ) + + m.fs.cv.add_momentum_balances( + balance_type=MomentumBalanceType.pressureTotal, has_pressure_change=True + ) + + m.fs.cv.apply_transformation() + + scaler_obj = ControlVolume1DScaler() + scaler_obj.default_scaling_factors["area"] = 1 / 83 + scaler_obj.default_scaling_factors["length"] = 1 / 599 + scaler_obj.scale_model(m.fs.cv) + + cv = m.fs.cv + + # Variables + assert get_scaling_factor(cv.length) == 1 / 599 + for v in cv.area.values(): + assert get_scaling_factor(v) == 1 / 83 + + for v in cv._enthalpy_flow.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.enthalpy_flow_dx.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.pressure_dx.values(): + assert get_scaling_factor(v) == 1 / 13 + + for v in cv.heat.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.work.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.enthalpy_transfer.values(): + assert get_scaling_factor(v) == approx(599 / 37) + + for v in cv.deltaP.values(): + assert get_scaling_factor(v) == approx(599 / 13) + + # Constraints + for v in cv.enthalpy_flow_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 37 + + for v in cv.pressure_dx_disc_eq.values(): + assert get_scaling_factor(v) == 1 / 13 + + for c in cv.enthalpy_balances.values(): + assert get_scaling_factor(c) == approx(1 / 37) + + for c in cv.pressure_balance.values(): + assert get_scaling_factor(c) == approx(1 / 13) + + assert len(list_unscaled_variables(m, include_fixed=True)) == 0 + assert len(list_unscaled_constraints(m)) == 0 + + +# TODO element balance +# @pytest.mark.unit +# def test_full_auto_scaling_mbtype_element(): +# m = pyo.ConcreteModel() +# m.fs = FlowsheetBlock(dynamic=True, time_units=pyo.units.s) +# m.fs.pp = PhysicalParameterTestBlock() +# m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp) +# m.fs.cv = ControlVolume1DBlock( +# property_package=m.fs.pp, +# reaction_package=m.fs.rp, +# transformation_method="dae.finite_difference", +# transformation_scheme="BACKWARD", +# finite_elements=10, +# ) +# m.fs.cv.add_geometry() +# m.fs.cv.add_state_blocks(has_phase_equilibrium=False) +# m.fs.cv.add_reaction_blocks(has_equilibrium=False) + +# m.fs.cv.add_total_element_balances(has_mass_transfer=True) + +# m.fs.cv.apply_transformation() +# m.discretizer = pyo.TransformationFactory("dae.finite_difference") +# m.discretizer.apply_to(m, nfe=3, wrt=m.fs.time, scheme="BACKWARD") + +# iscale.calculate_scaling_factors(m) + +# # check that all variables have scaling factors +# unscaled_var_list = list(iscale.unscaled_variables_generator(m)) +# # Unscaled variables are: +# # cp (44 space and time points) +# assert len(unscaled_var_list) == 44 +# # check that all constraints have been scaled +# unscaled_constraint_list = list(iscale.unscaled_constraints_generator(m)) +# assert len(unscaled_constraint_list) == 0 diff --git a/idaes/core/base/tests/test_control_volume_base.py b/idaes/core/base/tests/test_control_volume_base.py index 9bd2b53e6e..d39697d6e3 100644 --- a/idaes/core/base/tests/test_control_volume_base.py +++ b/idaes/core/base/tests/test_control_volume_base.py @@ -15,6 +15,7 @@ Author: Andrew Lee """ +import re import inspect import pytest from types import MethodType @@ -35,6 +36,7 @@ PhysicalParameterBlock, ReactionParameterBlock, ) +from idaes.core.base.control_volume_base import ControlVolumeScalerBase from idaes.core.util.exceptions import ( ConfigurationError, DynamicError, @@ -73,13 +75,36 @@ def test_momentum_balance_type(): @pytest.mark.unit def testflow_direction(): - assert len(FlowDirection) == 2 + assert len(FlowDirection) == 3 # Test that error is raised when given non-member with pytest.raises(AttributeError): FlowDirection.foo # pylint: disable=no-member +@pytest.mark.unit +def test_ControlVolumeScalerBase_no_state_block_ref(): + m = ConcreteModel() + scaler_obj = ControlVolumeScalerBase() + with pytest.raises( + AttributeError, + match=re.escape( + "The _state_block_ref attribute was not overridden by the " + "class inheriting from ControlVolumeScalerBase." + ), + ): + scaler_obj.scale_model(m) + + with pytest.raises( + AttributeError, + match=re.escape( + "The _state_block_ref attribute was not overridden by the " + "class inheriting from ControlVolumeScalerBase." + ), + ): + scaler_obj.constraint_scaling_routine(m) + + # ----------------------------------------------------------------------------- # Test CONFIG_Template @pytest.mark.unit diff --git a/idaes/core/scaling/custom_scaler_base.py b/idaes/core/scaling/custom_scaler_base.py index 49d35dd6c6..8aa7b0746a 100644 --- a/idaes/core/scaling/custom_scaler_base.py +++ b/idaes/core/scaling/custom_scaler_base.py @@ -299,6 +299,7 @@ def scale_variable_by_component( variable=target_variable, scaling_factor=sf, overwrite=overwrite ) else: + # TODO add infrastructure to log a warning _log.debug( f"Could not set scaling factor for {target_variable.name}, " f"no scaling factor set for {scaling_component.name}" @@ -893,8 +894,7 @@ def call_submodel_scaler_method( if submodel_scalers is None: submodel_scalers = {} - # Iterate over indices of submodel - for smdata in submodel.values(): + def scale_smdata(smdata): # Get Scaler for submodel if submodel in submodel_scalers: scaler = submodel_scalers[submodel] @@ -929,3 +929,10 @@ def call_submodel_scaler_method( f"Scaler for {submodel} does not have a method named {method}." ) from err smeth(smdata, submodel_scalers=submodel_scalers, overwrite=overwrite) + + if submodel.is_indexed(): + # Iterate over indices of submodel + for data in submodel.values(): + scale_smdata(data) + else: + scale_smdata(submodel) diff --git a/idaes/core/scaling/tests/test_custom_scaler_base.py b/idaes/core/scaling/tests/test_custom_scaler_base.py index 2bc4a20d02..1b7ccbc147 100644 --- a/idaes/core/scaling/tests/test_custom_scaler_base.py +++ b/idaes/core/scaling/tests/test_custom_scaler_base.py @@ -1117,6 +1117,71 @@ def test_propagate_state_scaling_scalar_to_indexed(self): assert sd.scaling_factor[j] == 20 * 2 * count count += 1 + @pytest.mark.unit + def test_propagate_state_scaling_scalar_to_indexed_within_block(self): + # Dummy up two state blocks + m = ConcreteModel() + + m.properties = PhysicalParameterTestBlock() + + m.state1 = m.properties.build_state_block([1, 2, 3]) + + # Set scaling factors on state1 + for t, sd in m.state1.items(): + sd.scaling_factor = Suffix(direction=Suffix.EXPORT) + sd.scaling_factor[sd.temperature] = 100 * t + sd.scaling_factor[sd.pressure] = 1e5 * t + + count = 1 + for j in sd.flow_mol_phase_comp.values(): + sd.scaling_factor[j] = 10 * t * count + count += 1 + + sb = CustomScalerBase() + sb.propagate_state_scaling(m.state1, m.state1[1]) + + for t, sd in m.state1.items(): + assert sd.scaling_factor[sd.temperature] == 100 * t + assert sd.scaling_factor[sd.pressure] == 1e5 * t + + count = 1 + for j in sd.flow_mol_phase_comp.values(): + assert sd.scaling_factor[j] == 10 * t * count + count += 1 + + # Test for overwrite=False + for t, sd in m.state1.items(): + sd.scaling_factor[sd.temperature] = 200 * t + sd.scaling_factor[sd.pressure] = 2e5 * t + + count = 1 + for j in sd.flow_mol_phase_comp.values(): + sd.scaling_factor[j] = 20 * t * count + count += 1 + + sb.propagate_state_scaling(m.state1, m.state1[1], overwrite=False) + + for t, sd in m.state1.items(): + assert sd.scaling_factor[sd.temperature] == 200 * t + assert sd.scaling_factor[sd.pressure] == 2e5 * t + + count = 1 + for j in sd.flow_mol_phase_comp.values(): + assert sd.scaling_factor[j] == 20 * t * count + count += 1 + + # Test for overwrite=True + sb.propagate_state_scaling(m.state1, m.state1[2], overwrite=True) + + for t, sd in m.state1.items(): + assert sd.scaling_factor[sd.temperature] == 200 * 2 + assert sd.scaling_factor[sd.pressure] == 2e5 * 2 + + count = 1 + for j in sd.flow_mol_phase_comp.values(): + assert sd.scaling_factor[j] == 20 * 2 * count + count += 1 + @pytest.mark.unit def test_propagate_state_scaling_scalar_to_scalar(self): # Dummy up two state blocks @@ -1266,5 +1331,26 @@ def test_call_submodel_scaler_method_user_scaler_class(self, caplog): assert "Using user-defined Scaler for b." in caplog.text + @pytest.mark.unit + def test_call_submodel_scaler_method_default_scaler_blockdata(self, caplog): + caplog.set_level(idaeslog.DEBUG, logger="idaes") + + # Dummy up a nested model + m = ConcreteModel() + m.b = Block([1, 2, 3]) + for bd in m.b.values(): + bd.default_scaler = DummyScaler + + sb = CustomScalerBase() + sb.call_submodel_scaler_method(m.b[1], method="dummy_method", overwrite=True) + + for i, bd in m.b.items(): + if i == 1: + assert bd._dummy_scaler_test + else: + assert not hasattr(bd, "_dummy_scaler_test") + + assert "Using default Scaler for b[1]." in caplog.text + # TODO additional tests for nested submodel scalers. diff --git a/idaes/models/properties/examples/saponification_thermo.py b/idaes/models/properties/examples/saponification_thermo.py index 7086858fe7..20f0ddf07d 100644 --- a/idaes/models/properties/examples/saponification_thermo.py +++ b/idaes/models/properties/examples/saponification_thermo.py @@ -60,7 +60,8 @@ class PhysicalParameterData(PhysicalParameterBlock): Property Parameter Block Class Contains parameters and indexing sets associated with properties for - superheated steam. + a dilute solution of NaOH, Ethyl Acetate, Sodium Acetate, and Ethanol + in water. """ diff --git a/idaes/models/unit_models/mixer.py b/idaes/models/unit_models/mixer.py index 6cd057963b..f9dd644856 100644 --- a/idaes/models/unit_models/mixer.py +++ b/idaes/models/unit_models/mixer.py @@ -18,6 +18,7 @@ from pyomo.environ import ( Block, check_optimal_termination, + ComponentMap, Constraint, Param, PositiveReals, @@ -34,6 +35,7 @@ MaterialBalanceType, MaterialFlowBasis, ) +from idaes.core.base.control_volume_base import ControlVolumeScalerBase from idaes.core.util.config import ( is_physical_parameter_block, is_state_block, @@ -52,7 +54,7 @@ import idaes.logger as idaeslog -__author__ = "Andrew Lee" +__author__ = "Andrew Lee, Douglas Allan" # Set up logger @@ -80,6 +82,96 @@ class MomentumMixingType(Enum): minimize_and_equality = 3 +class MixerScaler(ControlVolumeScalerBase): + """ + Scaler object for the Mixer unit model + """ + + # This attribute gives the parent ControlVolumeScalerBase + # methods a state block with the same index as the material + # and energy balances to get scaling information from + _state_block_ref = "mixed_state" + + def variable_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None + ): + for _, iblock in model.inlet_blocks: + self.call_submodel_scaler_method( + submodel=iblock, + submodel_scalers=submodel_scalers, + method="variable_scaling_routine", + overwrite=overwrite, + ) + + if model.config.mixed_state_block is None: + mblock = model.mixed_state + else: + mblock = model.config.mixed_state_block + + self.call_submodel_scaler_method( + submodel=mblock, + submodel_scalers=submodel_scalers, + method="variable_scaling_routine", + overwrite=overwrite, + ) + + # Scale inherent reaction and phase equilibrium + # variables, if they exist + super().variable_scaling_routine( + model, overwrite=overwrite, submodel_scalers=submodel_scalers + ) + + if hasattr(model, "minimum_pressure"): + for (t, _), v in model.minimum_pressure.items(): + self.scale_variable_by_component( + v, model.mixed_state[t].pressure, overwrite=overwrite + ) + + def constraint_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None + ): + for _, iblock in model.inlet_blocks: + self.call_submodel_scaler_method( + submodel=iblock, + submodel_scalers=submodel_scalers, + method="constraint_scaling_routine", + overwrite=overwrite, + ) + + if model.config.mixed_state_block is None: + mblock = model.mixed_state + else: + mblock = model.config.mixed_state_block + + self.call_submodel_scaler_method( + submodel=mblock, + submodel_scalers=submodel_scalers, + method="constraint_scaling_routine", + overwrite=overwrite, + ) + + # Scale material and energy balance equations + super().constraint_scaling_routine( + model, overwrite=overwrite, submodel_scalers=submodel_scalers + ) + + if hasattr(model, "pressure_equality_constraints"): + for (t, _), con in model.pressure_equality_constraints.items(): + self.scale_constraint_by_component( + con, model.mixed_state[t].pressure, overwrite=overwrite + ) + if hasattr(model, "minimum_pressure_constraint"): + for (t, _), con in model.minimum_pressure_constraint.items(): + self.scale_constraint_by_component( + con, model.mixed_state[t].pressure, overwrite=overwrite + ) + if hasattr(model, "mixture_pressure"): + for t, con in model.mixture_pressure.items(): + self.scale_constraint_by_component( + con, model.mixed_state[t].pressure, overwrite=overwrite + ) + + class MixerInitializer(ModularInitializerBase): """ Hierarchical Initializer for Mixer blocks. @@ -113,15 +205,11 @@ def initialization_routine( solver = self._get_solver() # Initialize inlet state blocks - inlet_list = model.create_inlet_list() i_block_list = [] - for i in inlet_list: - i_block = getattr(model, i + "_state") - i_block_list.append(i_block) - # Get initializer for inlet + for _, i_block in model.inlet_blocks: + i_block_list.append(i_block) iinit = self.get_submodel_initializer(i_block) - iinit.initialize(i_block) # Initialize mixed state block @@ -225,6 +313,9 @@ class MixerData(UnitModelBlockData): """ default_initializer = MixerInitializer + default_scaler = MixerScaler + + _inlet_dict = None CONFIG = ConfigBlock() CONFIG.declare( @@ -397,6 +488,19 @@ class MixerData(UnitModelBlockData): ), ) + @property + def inlet_blocks(self): + """ + Allows the user to iterate over the inlet stream names and state blocks. + + Returns: + dict_items with the inlet stream names as keys and the + inlet state blocks as values + """ + # Return an iterator so the user cannot + # mutate _inlet_dict + return self._inlet_dict.items() + def build(self): """ General build method for MixerData. This method calls a number @@ -424,6 +528,8 @@ def build(self): # Build StateBlocks inlet_blocks = self.add_inlet_state_blocks(inlet_list) + self._inlet_dict = {key: blk for key, blk in zip(inlet_list, inlet_blocks)} + if self.config.mixed_state_block is None: mixed_block = self.add_mixed_state_block() else: @@ -635,6 +741,7 @@ def add_material_mixing_equations(self, inlet_blocks, mixed_block, mb_type): else: # Let this pass for now with no units flow_units = None + self._constructed_material_balance_type = mb_type if mixed_block.include_inherent_reactions: if mb_type == MaterialBalanceType.total: diff --git a/idaes/models/unit_models/tests/test_mixer.py b/idaes/models/unit_models/tests/test_mixer.py index 1d15826efd..3494a4fd25 100644 --- a/idaes/models/unit_models/tests/test_mixer.py +++ b/idaes/models/unit_models/tests/test_mixer.py @@ -13,18 +13,20 @@ """ Tests for ControlVolumeBlockData. -Author: Andrew Lee +Author: Andrew Lee, Douglas Allan """ import pytest import pandas from pyomo.environ import ( check_optimal_termination, + ComponentMap, ConcreteModel, Constraint, Param, RangeSet, Set, + TransformationFactory, Var, value, units as pyunits, @@ -46,6 +48,11 @@ Component, MaterialFlowBasis, ) +from idaes.core.scaling import CustomScalerBase +from idaes.core.util.scaling import ( + get_jacobian, + jacobian_cond, +) from idaes.models.properties.activity_coeff_models.BTX_activity_coeff_VLE import ( BTXParameterBlock, ) @@ -60,6 +67,7 @@ MixingType, MomentumMixingType, MixerInitializer, + MixerScaler, ) from idaes.core.util.exceptions import ( BurntToast, @@ -593,6 +601,10 @@ def test_build_default(self): m.fs.pp = PhysicalParameterTestBlock() m.fs.mix = Mixer(property_package=m.fs.pp) + inlet_dict = {key: val for key, val in m.fs.mix.inlet_blocks} + assert len(inlet_dict) == 2 + assert inlet_dict["inlet_1"] is m.fs.mix.inlet_1_state + assert inlet_dict["inlet_2"] is m.fs.mix.inlet_2_state assert isinstance(m.fs.mix.material_mixing_equations, Constraint) assert len(m.fs.mix.material_mixing_equations) == 4 @@ -605,8 +617,9 @@ def test_build_default(self): assert isinstance(m.fs.mix.minimum_pressure, Var) assert len(m.fs.mix.minimum_pressure) == 2 assert isinstance(m.fs.mix.eps_pressure, Param) + assert len(m.fs.mix.minimum_pressure_constraint) == 2 assert isinstance(m.fs.mix.minimum_pressure_constraint, Constraint) - assert len(m.fs.mix.minimum_pressure) == 2 + assert len(m.fs.mix.mixture_pressure) == 1 assert isinstance(m.fs.mix.mixture_pressure, Constraint) assert isinstance(m.fs.mix.inlet_1, Port) @@ -633,8 +646,9 @@ def test_build_phase_equilibrium(self): assert isinstance(m.fs.mix.minimum_pressure, Var) assert len(m.fs.mix.minimum_pressure) == 2 assert isinstance(m.fs.mix.eps_pressure, Param) + assert len(m.fs.mix.minimum_pressure_constraint) == 2 assert isinstance(m.fs.mix.minimum_pressure_constraint, Constraint) - assert len(m.fs.mix.minimum_pressure) == 2 + assert len(m.fs.mix.mixture_pressure) == 1 assert isinstance(m.fs.mix.mixture_pressure, Constraint) assert isinstance(m.fs.mix.inlet_1, Port) @@ -665,6 +679,41 @@ def test_build_phase_pressure_equality(self): assert isinstance(m.fs.mix.inlet_2, Port) assert isinstance(m.fs.mix.outlet, Port) + @pytest.mark.build + @pytest.mark.unit + def test_build_named_inlets(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + + m.fs.mix = Mixer(property_package=m.fs.pp, inlet_list=["foo", "bar", "biz"]) + inlet_dict = {key: val for key, val in m.fs.mix.inlet_blocks} + assert len(inlet_dict) == 3 + assert inlet_dict["foo"] is m.fs.mix.foo_state + assert inlet_dict["bar"] is m.fs.mix.bar_state + assert inlet_dict["biz"] is m.fs.mix.biz_state + + assert isinstance(m.fs.mix.material_mixing_equations, Constraint) + assert len(m.fs.mix.material_mixing_equations) == 4 + assert hasattr(m.fs.mix, "phase_equilibrium_idx_ref") is False + + assert isinstance(m.fs.mix.enthalpy_mixing_equations, Constraint) + assert len(m.fs.mix.enthalpy_mixing_equations) == 1 + + assert m.fs.mix.inlet_idx.ctype is RangeSet + assert isinstance(m.fs.mix.minimum_pressure, Var) + assert len(m.fs.mix.minimum_pressure) == 3 + assert isinstance(m.fs.mix.eps_pressure, Param) + assert len(m.fs.mix.minimum_pressure_constraint) == 3 + assert isinstance(m.fs.mix.minimum_pressure_constraint, Constraint) + assert len(m.fs.mix.mixture_pressure) == 1 + assert isinstance(m.fs.mix.mixture_pressure, Constraint) + + assert isinstance(m.fs.mix.foo, Port) + assert isinstance(m.fs.mix.bar, Port) + assert isinstance(m.fs.mix.biz, Port) + assert isinstance(m.fs.mix.outlet, Port) + # ------------------------------------------------------------------------- # Test models checks, initialize and release state methods @pytest.mark.unit @@ -1242,61 +1291,185 @@ def test_numerical_issues(self, iapws): dt.assert_no_numerical_warnings() -# ----------------------------------------------------------------------------- -class TestSaponification(object): - @pytest.fixture(scope="class") - def sapon(self): - m = ConcreteModel() - m.fs = FlowsheetBlock(dynamic=False) +def _make_sapon_model(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) - m.fs.properties = SaponificationParameterBlock() + m.fs.properties = SaponificationParameterBlock() + + m.fs.unit = Mixer(property_package=m.fs.properties) + + m.fs.unit.inlet_1.flow_vol[0].fix(1e-3) + m.fs.unit.inlet_1.temperature[0].fix(320) + m.fs.unit.inlet_1.pressure[0].fix(101325) + m.fs.unit.inlet_1.conc_mol_comp[0, "H2O"].fix(55388.0) + m.fs.unit.inlet_1.conc_mol_comp[0, "NaOH"].fix(100.0) + m.fs.unit.inlet_1.conc_mol_comp[0, "EthylAcetate"].fix(100.0) + m.fs.unit.inlet_1.conc_mol_comp[0, "SodiumAcetate"].fix(1e-8) + m.fs.unit.inlet_1.conc_mol_comp[0, "Ethanol"].fix(1e-8) + + m.fs.unit.inlet_2.flow_vol[0].fix(1e-3) + m.fs.unit.inlet_2.temperature[0].fix(300) + m.fs.unit.inlet_2.pressure[0].fix(101325) + m.fs.unit.inlet_2.conc_mol_comp[0, "H2O"].fix(55388.0) + m.fs.unit.inlet_2.conc_mol_comp[0, "NaOH"].fix(100.0) + m.fs.unit.inlet_2.conc_mol_comp[0, "EthylAcetate"].fix(100.0) + m.fs.unit.inlet_2.conc_mol_comp[0, "SodiumAcetate"].fix(1e-8) + m.fs.unit.inlet_2.conc_mol_comp[0, "Ethanol"].fix(1e-8) + + return m + + +_sapon_expected = pandas.DataFrame.from_dict( + { + "Units": { + "Volumetric Flowrate": getattr(pyunits.pint_registry, "m**3/second"), + "Molar Concentration H2O": getattr(pyunits.pint_registry, "mole/m**3"), + "Molar Concentration NaOH": getattr(pyunits.pint_registry, "mole/m**3"), + "Molar Concentration EthylAcetate": getattr( + pyunits.pint_registry, "mole/m**3" + ), + "Molar Concentration SodiumAcetate": getattr( + pyunits.pint_registry, "mole/m**3" + ), + "Molar Concentration Ethanol": getattr(pyunits.pint_registry, "mole/m**3"), + "Temperature": getattr(pyunits.pint_registry, "K"), + "Pressure": getattr(pyunits.pint_registry, "Pa"), + }, + "inlet_1": { + "Volumetric Flowrate": 1e-3, + "Molar Concentration H2O": 55388, + "Molar Concentration NaOH": 100.00, + "Molar Concentration EthylAcetate": 100.00, + "Molar Concentration SodiumAcetate": 0, + "Molar Concentration Ethanol": 0, + "Temperature": 320, + "Pressure": 1.0132e05, + }, + "inlet_2": { + "Volumetric Flowrate": 1e-3, + "Molar Concentration H2O": 55388, + "Molar Concentration NaOH": 100.00, + "Molar Concentration EthylAcetate": 100.00, + "Molar Concentration SodiumAcetate": 0, + "Molar Concentration Ethanol": 0, + "Temperature": 300, + "Pressure": 1.0132e05, + }, + "Outlet": { + "Volumetric Flowrate": 1.00, + "Molar Concentration H2O": 100.00, + "Molar Concentration NaOH": 100.0, + "Molar Concentration EthylAcetate": 100.00, + "Molar Concentration SodiumAcetate": 100.00, + "Molar Concentration Ethanol": 100.00, + "Temperature": 298.15, + "Pressure": 1.0132e05, + }, + } +) - m.fs.unit = Mixer(property_package=m.fs.properties) - m.fs.unit.inlet_1.flow_vol[0].fix(1e-3) - m.fs.unit.inlet_1.temperature[0].fix(320) - m.fs.unit.inlet_1.pressure[0].fix(101325) - m.fs.unit.inlet_1.conc_mol_comp[0, "H2O"].fix(55388.0) - m.fs.unit.inlet_1.conc_mol_comp[0, "NaOH"].fix(100.0) - m.fs.unit.inlet_1.conc_mol_comp[0, "EthylAcetate"].fix(100.0) - m.fs.unit.inlet_1.conc_mol_comp[0, "SodiumAcetate"].fix(1e-8) - m.fs.unit.inlet_1.conc_mol_comp[0, "Ethanol"].fix(1e-8) - - m.fs.unit.inlet_2.flow_vol[0].fix(1e-3) - m.fs.unit.inlet_2.temperature[0].fix(300) - m.fs.unit.inlet_2.pressure[0].fix(101325) - m.fs.unit.inlet_2.conc_mol_comp[0, "H2O"].fix(55388.0) - m.fs.unit.inlet_2.conc_mol_comp[0, "NaOH"].fix(100.0) - m.fs.unit.inlet_2.conc_mol_comp[0, "EthylAcetate"].fix(100.0) - m.fs.unit.inlet_2.conc_mol_comp[0, "SodiumAcetate"].fix(1e-8) - m.fs.unit.inlet_2.conc_mol_comp[0, "Ethanol"].fix(1e-8) +def _test_build_sapon(sapon): + assert len(sapon.fs.unit.inlet_1.vars) == 4 + assert hasattr(sapon.fs.unit.inlet_1, "flow_vol") + assert hasattr(sapon.fs.unit.inlet_1, "conc_mol_comp") + assert hasattr(sapon.fs.unit.inlet_1, "temperature") + assert hasattr(sapon.fs.unit.inlet_1, "pressure") - return m + assert len(sapon.fs.unit.inlet_2.vars) == 4 + assert hasattr(sapon.fs.unit.inlet_2, "flow_vol") + assert hasattr(sapon.fs.unit.inlet_2, "conc_mol_comp") + assert hasattr(sapon.fs.unit.inlet_2, "temperature") + assert hasattr(sapon.fs.unit.inlet_2, "pressure") + + assert len(sapon.fs.unit.outlet.vars) == 4 + assert hasattr(sapon.fs.unit.outlet, "flow_vol") + assert hasattr(sapon.fs.unit.outlet, "conc_mol_comp") + assert hasattr(sapon.fs.unit.outlet, "temperature") + assert hasattr(sapon.fs.unit.outlet, "pressure") + + assert number_variables(sapon) == 26 + assert number_total_constraints(sapon) == 10 + assert number_unused_variables(sapon) == 0 + + +def _test_solution_sapon(sapon): + assert pytest.approx(2e-3, abs=1e-6) == value(sapon.fs.unit.outlet.flow_vol[0]) + + assert pytest.approx(55388.0, abs=1e0) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "H2O"] + ) + assert pytest.approx(100.0, abs=1e-3) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "NaOH"] + ) + assert pytest.approx(100.0, abs=1e-3) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "EthylAcetate"] + ) + assert pytest.approx(0.0, abs=1e-3) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "SodiumAcetate"] + ) + assert pytest.approx(0.0, abs=1e-3) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "Ethanol"] + ) + + assert pytest.approx(310.0, abs=1e-1) == value(sapon.fs.unit.outlet.temperature[0]) + + assert pytest.approx(101325, abs=1e2) == value(sapon.fs.unit.outlet.pressure[0]) + + +def _test_conservation_sapon(sapon): + assert ( + abs( + value( + sapon.fs.unit.inlet_1.flow_vol[0] + + sapon.fs.unit.inlet_2.flow_vol[0] + - sapon.fs.unit.outlet.flow_vol[0] + ) + ) + <= 1e-6 + ) + + assert ( + abs( + value( + sapon.fs.unit.inlet_1.flow_vol[0] + * sapon.fs.properties.dens_mol + * sapon.fs.properties.cp_mol + * ( + sapon.fs.unit.inlet_1.temperature[0] + - sapon.fs.properties.temperature_ref + ) + + sapon.fs.unit.inlet_2.flow_vol[0] + * sapon.fs.properties.dens_mol + * sapon.fs.properties.cp_mol + * ( + sapon.fs.unit.inlet_2.temperature[0] + - sapon.fs.properties.temperature_ref + ) + - sapon.fs.unit.outlet.flow_vol[0] + * sapon.fs.properties.dens_mol + * sapon.fs.properties.cp_mol + * ( + sapon.fs.unit.outlet.temperature[0] + - sapon.fs.properties.temperature_ref + ) + ) + ) + <= 1e-3 + ) + + +# ----------------------------------------------------------------------------- +class TestSaponificationLegacyScaling(object): + @pytest.fixture(scope="class") + def sapon(self): + return _make_sapon_model() @pytest.mark.build @pytest.mark.unit def test_build(self, sapon): - assert len(sapon.fs.unit.inlet_1.vars) == 4 - assert hasattr(sapon.fs.unit.inlet_1, "flow_vol") - assert hasattr(sapon.fs.unit.inlet_1, "conc_mol_comp") - assert hasattr(sapon.fs.unit.inlet_1, "temperature") - assert hasattr(sapon.fs.unit.inlet_1, "pressure") - - assert len(sapon.fs.unit.inlet_2.vars) == 4 - assert hasattr(sapon.fs.unit.inlet_2, "flow_vol") - assert hasattr(sapon.fs.unit.inlet_2, "conc_mol_comp") - assert hasattr(sapon.fs.unit.inlet_2, "temperature") - assert hasattr(sapon.fs.unit.inlet_2, "pressure") - - assert len(sapon.fs.unit.outlet.vars) == 4 - assert hasattr(sapon.fs.unit.outlet, "flow_vol") - assert hasattr(sapon.fs.unit.outlet, "conc_mol_comp") - assert hasattr(sapon.fs.unit.outlet, "temperature") - assert hasattr(sapon.fs.unit.outlet, "pressure") - - assert number_variables(sapon) == 26 - assert number_total_constraints(sapon) == 10 - assert number_unused_variables(sapon) == 0 + _test_build_sapon(sapon) @pytest.mark.component def test_structural_issues(self, sapon): @@ -1308,64 +1481,7 @@ def test_structural_issues(self, sapon): def test_get_stream_table_contents(self, sapon): stable = sapon.fs.unit._get_stream_table_contents() - expected = pandas.DataFrame.from_dict( - { - "Units": { - "Volumetric Flowrate": getattr( - pyunits.pint_registry, "m**3/second" - ), - "Molar Concentration H2O": getattr( - pyunits.pint_registry, "mole/m**3" - ), - "Molar Concentration NaOH": getattr( - pyunits.pint_registry, "mole/m**3" - ), - "Molar Concentration EthylAcetate": getattr( - pyunits.pint_registry, "mole/m**3" - ), - "Molar Concentration SodiumAcetate": getattr( - pyunits.pint_registry, "mole/m**3" - ), - "Molar Concentration Ethanol": getattr( - pyunits.pint_registry, "mole/m**3" - ), - "Temperature": getattr(pyunits.pint_registry, "K"), - "Pressure": getattr(pyunits.pint_registry, "Pa"), - }, - "inlet_1": { - "Volumetric Flowrate": 1e-3, - "Molar Concentration H2O": 55388, - "Molar Concentration NaOH": 100.00, - "Molar Concentration EthylAcetate": 100.00, - "Molar Concentration SodiumAcetate": 0, - "Molar Concentration Ethanol": 0, - "Temperature": 320, - "Pressure": 1.0132e05, - }, - "inlet_2": { - "Volumetric Flowrate": 1e-3, - "Molar Concentration H2O": 55388, - "Molar Concentration NaOH": 100.00, - "Molar Concentration EthylAcetate": 100.00, - "Molar Concentration SodiumAcetate": 0, - "Molar Concentration Ethanol": 0, - "Temperature": 300, - "Pressure": 1.0132e05, - }, - "Outlet": { - "Volumetric Flowrate": 1.00, - "Molar Concentration H2O": 100.00, - "Molar Concentration NaOH": 100.0, - "Molar Concentration EthylAcetate": 100.00, - "Molar Concentration SodiumAcetate": 100.00, - "Molar Concentration Ethanol": 100.00, - "Temperature": 298.15, - "Pressure": 1.0132e05, - }, - } - ) - - pandas.testing.assert_frame_equal(stable, expected, rtol=1e-4, atol=1e-4) + pandas.testing.assert_frame_equal(stable, _sapon_expected, rtol=1e-4, atol=1e-4) @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @@ -1389,74 +1505,101 @@ def test_solve(self, sapon): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_solution(self, sapon): - assert pytest.approx(2e-3, abs=1e-6) == value(sapon.fs.unit.outlet.flow_vol[0]) + _test_solution_sapon(sapon) - assert pytest.approx(55388.0, abs=1e0) == value( - sapon.fs.unit.outlet.conc_mol_comp[0, "H2O"] - ) - assert pytest.approx(100.0, abs=1e-3) == value( - sapon.fs.unit.outlet.conc_mol_comp[0, "NaOH"] - ) - assert pytest.approx(100.0, abs=1e-3) == value( - sapon.fs.unit.outlet.conc_mol_comp[0, "EthylAcetate"] - ) - assert pytest.approx(0.0, abs=1e-3) == value( - sapon.fs.unit.outlet.conc_mol_comp[0, "SodiumAcetate"] - ) - assert pytest.approx(0.0, abs=1e-3) == value( - sapon.fs.unit.outlet.conc_mol_comp[0, "Ethanol"] - ) + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_conservation(self, sapon): + _test_conservation_sapon(sapon) - assert pytest.approx(310.0, abs=1e-1) == value( - sapon.fs.unit.outlet.temperature[0] - ) + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_numerical_issues(self, sapon): + dt = DiagnosticsToolbox(sapon) + dt.assert_no_numerical_warnings() + + +class TestSaponificationScalerObject(object): + @pytest.fixture(scope="class") + def sapon(self): + return _make_sapon_model() - assert pytest.approx(101325, abs=1e2) == value(sapon.fs.unit.outlet.pressure[0]) + @pytest.mark.build + @pytest.mark.unit + def test_build(self, sapon): + _test_build_sapon(sapon) + + @pytest.mark.component + def test_structural_issues(self, sapon): + dt = DiagnosticsToolbox(sapon) + dt.assert_no_structural_warnings() + + @pytest.mark.ui + @pytest.mark.unit + def test_get_stream_table_contents(self, sapon): + stable = sapon.fs.unit._get_stream_table_contents() + + pandas.testing.assert_frame_equal(stable, _sapon_expected, rtol=1e-4, atol=1e-4) @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component - def test_conservation(self, sapon): - assert ( - abs( - value( - sapon.fs.unit.inlet_1.flow_vol[0] - + sapon.fs.unit.inlet_2.flow_vol[0] - - sapon.fs.unit.outlet.flow_vol[0] - ) - ) - <= 1e-6 - ) + def test_initialize(self, sapon): + initialization_tester(sapon) - assert ( - abs( - value( - sapon.fs.unit.inlet_1.flow_vol[0] - * sapon.fs.properties.dens_mol - * sapon.fs.properties.cp_mol - * ( - sapon.fs.unit.inlet_1.temperature[0] - - sapon.fs.properties.temperature_ref - ) - + sapon.fs.unit.inlet_2.flow_vol[0] - * sapon.fs.properties.dens_mol - * sapon.fs.properties.cp_mol - * ( - sapon.fs.unit.inlet_2.temperature[0] - - sapon.fs.properties.temperature_ref - ) - - sapon.fs.unit.outlet.flow_vol[0] - * sapon.fs.properties.dens_mol - * sapon.fs.properties.cp_mol - * ( - sapon.fs.unit.outlet.temperature[0] - - sapon.fs.properties.temperature_ref - ) - ) - ) - <= 1e-3 + @pytest.mark.component + def test_scaling(self, sapon): + unit = sapon.fs.unit + assert unit.default_scaler is MixerScaler + scaler_obj = MixerScaler() + scaler_obj.scale_model(unit) + + assert len(unit.scaling_factor) == 11 + + # Variables + for v in unit.minimum_pressure.values(): + assert unit.scaling_factor[v] == 1e-5 + + # Constraints + for idx, c in unit.material_mixing_equations.items(): + if idx[-1] == "H2O": + assert unit.scaling_factor[c] == 1e-2 + else: + assert unit.scaling_factor[c] == 1 + + assert unit.scaling_factor[unit.enthalpy_mixing_equations[0]] == pytest.approx( + 1.917e-6, rel=1e-3 ) + for c in unit.minimum_pressure_constraint.values(): + assert unit.scaling_factor[c] == 1e-5 + + assert unit.scaling_factor[unit.mixture_pressure[0]] == 1e-5 + # Expressions + assert not hasattr(unit, "scaling_hint") + + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solve(self, sapon): + results = solver.solve(sapon) + + # Check for optimal solution + assert check_optimal_termination(results) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solution(self, sapon): + _test_solution_sapon(sapon) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_conservation(self, sapon): + _test_conservation_sapon(sapon) + @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component @@ -1742,10 +1885,56 @@ def fix_initialization_states(self): self[k].hso4_dissociation.deactivate() +class LeachSolutionScaler(CustomScalerBase): + """ + Scaler object for leach solution properties + """ + + DEFAULT_SCALING_FACTORS = { + "flow_vol": 0.1, + "conc_mass_comp[H2O]": 1e-6, + "conc_mass_comp[H]": 1e-1, + "conc_mass_comp[HSO4]": 1e-2, + "conc_mass_comp[SO4]": 1e-1, + } + + def variable_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None + ): + self.scale_variable_by_default(model.flow_vol, overwrite=overwrite) + for var in model.conc_mass_comp.values(): + self.scale_variable_by_default(var, overwrite=overwrite) + + for idx, var in model.conc_mol_comp.items(): + self.scale_variable_by_definition_constraint( + var, model.molar_concentration_constraint[idx], overwrite=overwrite + ) + + def constraint_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: ComponentMap = None + ): + for idx, con in model.molar_concentration_constraint.items(): + # Constraint is on a mass basis + self.scale_constraint_by_component( + con, model.conc_mass_comp[idx], overwrite=overwrite + ) + if not model.config.defined_state: + self.scale_constraint_by_component( + model.h2o_concentration, + model.conc_mass_comp["H2O"], + overwrite=overwrite, + ) + self.scale_constraint_by_nominal_value( + model.hso4_dissociation, overwrite=overwrite + ) + + @declare_process_block_class( "LeachSolutionStateBlock", block_class=_LeachSolutionStateBlock ) class LeachSolutionStateBlockData(StateBlockData): + default_scaler = LeachSolutionScaler + def build(self): super().build() @@ -1809,62 +1998,153 @@ def define_state_vars(self): } -@pytest.mark.component -@pytest.mark.solver -def test_component_phase_w_inherent_rxns(): - m = ConcreteModel() - m.fs = FlowsheetBlock(dynamic=False) +class TestComponentPhaseWithInherentReactions(object): + @pytest.fixture(scope="class") + def model(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) - m.fs.properties = LeachSolutionParameters() + m.fs.properties = LeachSolutionParameters() - m.fs.unit = Mixer( - property_package=m.fs.properties, - material_balance_type=MaterialBalanceType.componentPhase, - energy_mixing_type=MixingType.none, - momentum_mixing_type=MomentumMixingType.none, - ) + m.fs.unit = Mixer( + property_package=m.fs.properties, + material_balance_type=MaterialBalanceType.componentPhase, + energy_mixing_type=MixingType.none, + momentum_mixing_type=MomentumMixingType.none, + ) - m.fs.unit.inlet_1.flow_vol[0].fix(10) - m.fs.unit.inlet_1.conc_mass_comp[0, "H2O"].fix(1e6 * pyunits.mg / pyunits.L) - m.fs.unit.inlet_1.conc_mass_comp[0, "H"].fix(57.54848 * pyunits.mg / pyunits.L) - m.fs.unit.inlet_1.conc_mass_comp[0, "HSO4"].fix(4117.798 * pyunits.mg / pyunits.L) - m.fs.unit.inlet_1.conc_mass_comp[0, "SO4"].fix(724.6539 * pyunits.mg / pyunits.L) + m.fs.unit.inlet_1.flow_vol[0].fix(10) + m.fs.unit.inlet_1.conc_mass_comp[0, "H2O"].fix(1e6 * pyunits.mg / pyunits.L) + m.fs.unit.inlet_1.conc_mass_comp[0, "H"].fix(57.54848 * pyunits.mg / pyunits.L) + m.fs.unit.inlet_1.conc_mass_comp[0, "HSO4"].fix( + 4117.798 * pyunits.mg / pyunits.L + ) + m.fs.unit.inlet_1.conc_mass_comp[0, "SO4"].fix( + 724.6539 * pyunits.mg / pyunits.L + ) - m.fs.unit.inlet_2.flow_vol[0].fix(10) - m.fs.unit.inlet_2.conc_mass_comp[0, "H2O"].fix(1e6 * pyunits.mg / pyunits.L) - m.fs.unit.inlet_2.conc_mass_comp[0, "H"].fix(57.54848 * pyunits.mg / pyunits.L) - m.fs.unit.inlet_2.conc_mass_comp[0, "HSO4"].fix(4117.798 * pyunits.mg / pyunits.L) - m.fs.unit.inlet_2.conc_mass_comp[0, "SO4"].fix(724.6539 * pyunits.mg / pyunits.L) + m.fs.unit.inlet_2.flow_vol[0].fix(10) + m.fs.unit.inlet_2.conc_mass_comp[0, "H2O"].fix(1e6 * pyunits.mg / pyunits.L) + m.fs.unit.inlet_2.conc_mass_comp[0, "H"].fix(57.54848 * pyunits.mg / pyunits.L) + m.fs.unit.inlet_2.conc_mass_comp[0, "HSO4"].fix( + 4117.798 * pyunits.mg / pyunits.L + ) + m.fs.unit.inlet_2.conc_mass_comp[0, "SO4"].fix( + 724.6539 * pyunits.mg / pyunits.L + ) - assert degrees_of_freedom(m) == 0 - assert_units_consistent(m) + return m - initializer = BlockTriangularizationInitializer() - initializer.initialize(m.fs.unit) + @pytest.mark.component + def test_structural_warnings(self, model): + dt = DiagnosticsToolbox(model) + dt.assert_no_structural_warnings() - results = solver.solve(m) - assert check_optimal_termination(results) + @pytest.mark.component + def test_scaling(self, model): + unit = model.fs.unit + assert unit.default_scaler is MixerScaler + MixerScaler().scale_model(unit) + + assert len(unit.scaling_factor) == 13 + + # Variables + assert unit.scaling_factor[ + unit.inherent_reaction_extent[0.0, "Ka2"] + ] == pytest.approx(960) + assert unit.scaling_factor[ + unit.inherent_reaction_generation[0.0, "liquid", "H2O"] + ] == pytest.approx(1.8e-3) + assert unit.scaling_factor[ + unit.inherent_reaction_generation[0.0, "liquid", "H"] + ] == pytest.approx(10) + assert unit.scaling_factor[ + unit.inherent_reaction_generation[0.0, "liquid", "HSO4"] + ] == pytest.approx(97) + assert unit.scaling_factor[ + unit.inherent_reaction_generation[0.0, "liquid", "SO4"] + ] == pytest.approx(960) + + # Constraints + assert unit.scaling_factor[ + unit.inherent_reaction_constraint[0.0, "liquid", "H2O"] + ] == pytest.approx(1.8e-3) + assert unit.scaling_factor[ + unit.inherent_reaction_constraint[0.0, "liquid", "H"] + ] == pytest.approx(10) + assert unit.scaling_factor[ + unit.inherent_reaction_constraint[0.0, "liquid", "HSO4"] + ] == pytest.approx(97) + assert unit.scaling_factor[ + unit.inherent_reaction_constraint[0.0, "liquid", "SO4"] + ] == pytest.approx(960) + assert unit.scaling_factor[ + unit.material_mixing_equations[0.0, "liquid", "H2O"] + ] == pytest.approx(1.8e-3) + assert unit.scaling_factor[ + unit.material_mixing_equations[0.0, "liquid", "H"] + ] == pytest.approx(10) + assert unit.scaling_factor[ + unit.material_mixing_equations[0.0, "liquid", "HSO4"] + ] == pytest.approx(97) + assert unit.scaling_factor[ + unit.material_mixing_equations[0.0, "liquid", "SO4"] + ] == pytest.approx(960) + + # Expressions + assert not hasattr(unit, "scaling_hint") - assert initializer.summary[m.fs.unit]["status"] == InitializationStatus.Ok + @pytest.mark.component + @pytest.mark.solver + def test_initialization(self, model): + initializer = BlockTriangularizationInitializer() + initializer.initialize(model.fs.unit) + assert initializer.summary[model.fs.unit]["status"] == InitializationStatus.Ok - assert value(m.fs.unit.outlet.flow_vol[0]) == pytest.approx(20, rel=1e-6) + @pytest.mark.component + @pytest.mark.solver + def test_solve(self, model): + results = solver.solve(model) + assert check_optimal_termination(results) - assert value(m.fs.unit.outlet.conc_mass_comp[0, "H2O"]) == pytest.approx( - 1e6, rel=1e-6 - ) - assert value(m.fs.unit.outlet.conc_mass_comp[0, "H"]) == pytest.approx( - 57.54848, rel=1e-6 - ) - assert value(m.fs.unit.outlet.conc_mass_comp[0, "HSO4"]) == pytest.approx( - 4117.798, rel=1e-6 - ) - assert value(m.fs.unit.outlet.conc_mass_comp[0, "SO4"]) == pytest.approx( - 724.6539, rel=1e-6 - ) + # Mark this test with solver because it relies on the + # results of the previous test, which is also marked + # solver. + @pytest.mark.component + @pytest.mark.solver + def test_solution(self, model): + unit = model.fs.unit + assert value(unit.outlet.flow_vol[0]) == pytest.approx(20, rel=1e-6) - assert value(m.fs.unit.inherent_reaction_extent[0, "Ka2"]) == pytest.approx( - 0, abs=1e-6 - ) + assert value(unit.outlet.conc_mass_comp[0, "H2O"]) == pytest.approx( + 1e6, rel=1e-6 + ) + assert value(unit.outlet.conc_mass_comp[0, "H"]) == pytest.approx( + 57.54848, rel=1e-6 + ) + assert value(unit.outlet.conc_mass_comp[0, "HSO4"]) == pytest.approx( + 4117.798, rel=1e-6 + ) + assert value(unit.outlet.conc_mass_comp[0, "SO4"]) == pytest.approx( + 724.6539, rel=1e-6 + ) + + assert value(unit.inherent_reaction_extent[0, "Ka2"]) == pytest.approx( + 0, abs=1e-6 + ) + + @pytest.mark.component + @pytest.mark.solver + def test_no_numerical_warnings(self, model): + # TODO revisit once the diagnostics toolbox takes scaling into account + sm = TransformationFactory("core.scale_model").create_using(model, rename=False) + dt = DiagnosticsToolbox(sm) + dt.assert_no_numerical_warnings() + + jac, _ = get_jacobian(sm, scaled=False) + assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx( + 1.074e4, rel=1e-3 + ) @pytest.mark.component