Mea column improvements

idaes/core/util/initialization.py

@@ -21,6 +21,7 @@
     Constraint,
     value,
 )
+from pyomo.common.collections import ComponentMap
 from pyomo.network import Arc
 from pyomo.dae import ContinuousSet
 from pyomo.core.expr.visitor import identify_variables
@@ -533,3 +534,35 @@ def initialize_by_time_element(fs, time, **kwargs):
 
     # Logger message that initialization is finished
     init_log.info("Initialization completed. Model has been reactivated")
+
+def _fix_vars(var_list):
+    flags = ComponentMap()
+    for var in var_list:
+        if var.is_indexed():
+            for subvar in var.values():
+                flags[subvar] = subvar.fixed
+                subvar.fix()
+        else:
+            flags[var] = var.fixed
+            var.fix()
+    return flags
+
+
+def _unfix_vars(var_list):
+    flags = ComponentMap()
+    if var.is_indexed():
+        for subvar in var.values():
+            flags[subvar] = subvar.fixed
+            subvar.unfix()
+    else:
+        flags[var] = var.fixed
+        var.unfix()
+    return flags
+
+
+def _restore_fixedness(flags):
+    for var, fix in flags.items():
+        if fix:
+            var.fix()
+        else:
+            var.unfix()

idaes/models/properties/modular_properties/base/generic_property.py

@@ -4689,7 +4689,7 @@ def _log_mole_frac_phase_comp(self):
         try:
             self.log_mole_frac_phase_comp = Var(
                 self.phase_component_set,
-                initialize=0,
+                initialize=-1,
                 bounds=(-100, log(1.001)),
                 units=pyunits.dimensionless,
                 doc="Log of mole fractions of component by phase",
@@ -4714,8 +4714,8 @@ def _log_mole_frac_phase_comp_apparent(self):
         try:
             self.log_mole_frac_phase_comp_apparent = Var(
                 self.params.apparent_phase_component_set,
-                initialize=0,
-                bounds=(-50, 0),
+                initialize=-1,
+                bounds=(-100, log(1.001)),
                 units=pyunits.dimensionless,
                 doc="Log of mole fractions of component by phase",
             )
@@ -4739,8 +4739,8 @@ def _log_mole_frac_phase_comp_true(self):
         try:
             self.log_mole_frac_phase_comp_true = Var(
                 self.params.true_phase_component_set,
-                initialize=0,
-                bounds=(-50, 0),
+                initialize=-1,
+                bounds=(-100, log(1.001)),
                 units=pyunits.dimensionless,
                 doc="Log of mole fractions of component by phase",
             )

idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py

@@ -633,8 +633,10 @@ def test_equilibrium_form(self, model):
         )
         assert value(model.rblock[1].equilibrium_constraint["e2"].body) == value(
             model.rblock[1].log_k_eq["e2"]
-            - model.sblock[1].log_mole_frac_phase_comp["p2", "c1"] * -5
-            + model.sblock[1].log_mole_frac_phase_comp["p2", "c2"] * 6
+            - (
+                model.sblock[1].log_mole_frac_phase_comp["p2", "c1"] * -5
+                + model.sblock[1].log_mole_frac_phase_comp["p2", "c2"] * 6
+            )
         )
 
     @pytest.mark.unit

idaes/models/properties/modular_properties/state_definitions/electrolyte_states.py

@@ -118,7 +118,7 @@ def appr_to_true_species(b, p, j):
             # Next, check for inherent reactions
             if b.has_inherent_reactions:
                 for r in b.params.inherent_reaction_idx:
-                    # Get stoichiometric coeffiicient for inherent reactions
+                    # Get stoichiometric coefficient for inherent reactions
                     gamma = b.params.inherent_reaction_stoichiometry[r, p, j]
 
                     if gamma != 0:
@@ -156,16 +156,37 @@ def true_species_mole_fractions(b, p, j):
 
 def _apparent_species_scaling(b):
     for p, j in b.params.true_phase_component_set:
-        sf = iscale.get_scaling_factor(
-            b.flow_mol_phase_comp_true[p, j], default=1, warning=True
-        )
+        pobj = b.params.get_phase(p)
+        cobj = b.params.get_component(j)
 
-        iscale.constraint_scaling_transform(
-            b.appr_to_true_species[p, j], sf, overwrite=False
-        )
-        iscale.constraint_scaling_transform(
-            b.true_mole_frac_constraint[p, j], sf, overwrite=False
-        )
+        # For apparent species in an aqueous phase, want to scale material balance
+        # equations by apparent magnitude. If they're not consumed to exhaustion,
+        # they'll be well-scaled in this equation. If they are consumed to "exhaustion"
+        # (several orders of magnitude below apparent concentration), then their concentration
+        # will decouple from this equation and will be given by equilibrium equation instead
+        if pobj.is_aqueous_phase and not isinstance(cobj, IonData):
+            sf_apparent = iscale.get_scaling_factor(
+                b.flow_mol_phase_comp_apparent[p, j], default=1, warning=True
+            )
+            iscale.constraint_scaling_transform(
+                b.appr_to_true_species[p, j], sf_apparent, overwrite=False
+            )
+            sf_true = iscale.get_scaling_factor(
+                b.flow_mol_phase_comp_true[p, j], default=1, warning=True
+            )
+            iscale.constraint_scaling_transform(
+                b.true_mole_frac_constraint[p, j], sf_true, overwrite=False
+            )
+        else:
+            sf = iscale.get_scaling_factor(
+                b.flow_mol_phase_comp_true[p, j], default=1, warning=True
+            )
+            iscale.constraint_scaling_transform(
+                b.appr_to_true_species[p, j], sf, overwrite=False
+            )
+            iscale.constraint_scaling_transform(
+                b.true_mole_frac_constraint[p, j], sf, overwrite=False
+            )
 
 
 def _true_species_state(b):

idaes/models/unit_models/heat_exchanger_ntu.py

@@ -48,6 +48,7 @@
 from idaes.core.util.math import smooth_min, smooth_max
 from idaes.core.solvers import get_solver
 from idaes.core.util.exceptions import InitializationError
+import idaes.core.util.scaling as iscale
 import idaes.logger as idaeslog
 from idaes.core.initialization import SingleControlVolumeUnitInitializer
 
@@ -478,7 +479,13 @@ def rule_entu(blk, t):
 
         self.heat_duty_constraint = Constraint(self.flowsheet().time, rule=rule_entu)
 
-    # TODO : Add scaling methods
+    def calculate_scaling_factors(self):
+        super().calculate_scaling_factors()
+
+        for t in self.flowsheet().time:
+            sf_heat = iscale.get_scaling_factor(self.hot_side.heat[t], default=1, warning=True)
+            iscale.constraint_scaling_transform(self.energy_balance_constraint[t], sf_heat, overwrite=False)
+            iscale.constraint_scaling_transform(self.heat_duty_constraint[t], sf_heat, overwrite=False)
 
     def initialize_build(
         self,