Scaling, fix minor MCAS issue, and set up Seawater properties for dynamics

watertap/core/membrane_channel0d.py

@@ -387,13 +387,18 @@ def calculate_scaling_factors(self):
                 if iscale.get_scaling_factor(v) is None:
                     iscale.set_scaling_factor(v, 1e3)
 
-        super().calculate_scaling_factors()
-
         if hasattr(self, "area"):
             if iscale.get_scaling_factor(self.area) is None:
                 iscale.set_scaling_factor(self.area, 100)
 
+        super().calculate_scaling_factors()
+
         if hasattr(self, "dP_dx"):
             for v in self.dP_dx.values():
                 if iscale.get_scaling_factor(v) is None:
                     iscale.set_scaling_factor(v, 1e-4)
+
+        if hasattr(self, "eq_pressure_change"):
+            for t, condata in self.eq_pressure_change.items():
+                sf = iscale.get_scaling_factor(self.deltaP[t], default=1)
+                iscale.constraint_scaling_transform(condata, sf)

watertap/core/membrane_channel1d.py

@@ -242,14 +242,14 @@ def initialize(
             self.release_state(source_flags, outlvl)
 
     def calculate_scaling_factors(self):
+        if iscale.get_scaling_factor(self.area) is None:
+            iscale.set_scaling_factor(self.area, 100)
+
         super().calculate_scaling_factors()
 
         # setting scaling factors for variables
 
         # will not override if the user provides the scaling factor
-        ## default of 1 set by ControlVolume1D
-        if iscale.get_scaling_factor(self.area) == 1:
-            iscale.set_scaling_factor(self.area, 100)
 
         if hasattr(self, "pressure_change_total"):
             for v in self.pressure_change_total.values():

watertap/core/membrane_channel_base.py

@@ -19,6 +19,7 @@
     Expression,
     Param,
     NonNegativeReals,
+    value,
     Var,
     units as pyunits,
 )
@@ -996,7 +997,11 @@ def calculate_scaling_factors(self):
         if hasattr(self, "velocity"):
             for v in self.velocity.values():
                 if iscale.get_scaling_factor(v) is None:
-                    iscale.set_scaling_factor(v, 1)
+                    sf_A = iscale.get_scaling_factor(self.area, default=1)
+                    sf_F_vol = iscale.get_scaling_factor(
+                        self.properties[t, x].flow_vol_phase["Liq"], default=1
+                    )
+                    iscale.set_scaling_factor(v, sf_F_vol / sf_A, 1)
 
         if hasattr(self, "friction_factor_darcy"):
             for v in self.friction_factor_darcy.values():
@@ -1008,6 +1013,71 @@ def calculate_scaling_factors(self):
                 if iscale.get_scaling_factor(v) is None:
                     iscale.set_scaling_factor(v, 1)
 
+        # TODO why did this file have zero scaling factors for constraints
+        if hasattr(self, "eq_dh"):
+            sf_dh = iscale.get_scaling_factor(self.dh, default=1)
+            iscale.constraint_scaling_transform(self.eq_dh, sf_dh)
+
+        if hasattr(self, "eq_area"):
+            sf_A = iscale.get_scaling_factor(self.area, default=1)
+            iscale.constraint_scaling_transform(self.eq_area, sf_A)
+
+        if hasattr(self, "eq_velocity"):
+            for (t, x), condata in self.eq_velocity.items():
+                sf_F_vol = iscale.get_scaling_factor(
+                    self.properties[t, x].flow_vol_phase["Liq"], default=1
+                )
+                iscale.constraint_scaling_transform(condata, sf_F_vol)
+
+        if hasattr(self, "eq_equal_flow_vol_interface"):
+            for (t, x), condata in self.eq_equal_flow_vol_interface.items():
+                sf_F_vol = min(
+                    iscale.get_scaling_factor(
+                        self.properties[t, x].flow_vol_phase["Liq"], default=1
+                    ),
+                    iscale.get_scaling_factor(
+                        self.properties_interface[t, x].flow_vol_phase["Liq"], default=1
+                    ),
+                )
+                iscale.constraint_scaling_transform(condata, sf_F_vol)
+
+        if hasattr(self, "eq_equal_pressure_interface"):
+            for (t, x), condata in self.eq_equal_pressure_interface.items():
+                sf_P = min(
+                    iscale.get_scaling_factor(
+                        self.properties_interface[t, x].pressure, default=1
+                    ),
+                    iscale.get_scaling_factor(
+                        self.properties[t, x].pressure, default=1
+                    ),
+                )
+                iscale.constraint_scaling_transform(condata, sf_P)
+
+        if hasattr(self, "eq_K"):
+            for (t, x, j), condata in self.eq_K.items():
+                sf_K = iscale.get_scaling_factor(self.K[t, x, j])
+                # Hopefully h has been calculated before this
+                sf_h = 1 / value(self.dh)
+                iscale.constraint_scaling_transform(condata, sf_K * sf_h)
+
+        if hasattr(self, "eq_N_Re"):
+            for (t, x), condata in self.eq_N_Re.items():
+                sf_Re = iscale.get_scaling_factor(self.N_Re[t, x], default=1)
+                sf_A = iscale.get_scaling_factor(
+                    self.area, default=1 / value(self.area)
+                )
+                sf_visc = iscale.get_scaling_factor(
+                    self.properties[t, x].visc_d_phase["Liq"], default=1
+                )
+                iscale.constraint_scaling_transform(condata, sf_Re * sf_A * sf_visc)
+
+        if hasattr(self, "eq_N_Sc_comp"):
+            for (t, x, j), condata in self.eq_N_Sc_comp.items():
+                sf_visc = iscale.get_scaling_factor(
+                    self.properties[t, x].visc_d_phase["Liq"], default=1
+                )
+                iscale.constraint_scaling_transform(condata, sf_visc)
+
 
 # helper for validating configuration arguments for this CV
 def validate_membrane_config_args(unit):

watertap/core/zero_order_sido.py

@@ -282,16 +282,17 @@ def initialize_sido(
 
 def calculate_scaling_factors_sido(self):
     # Get default scale factors and do calculations from base classes
-    for t, v in self.water_recovery_equation.items():
-        iscale.constraint_scaling_transform(
-            v,
-            iscale.get_scaling_factor(
-                self.properties_in[t].flow_mass_comp["H2O"],
-                default=1,
-                warning=True,
-                hint=" for water recovery",
-            ),
-        )
+    if hasattr(self, "water_recovery_equation"):
+        for t, v in self.water_recovery_equation.items():
+            iscale.constraint_scaling_transform(
+                v,
+                iscale.get_scaling_factor(
+                    self.properties_in[t].flow_mass_comp["H2O"],
+                    default=1,
+                    warning=True,
+                    hint=" for water recovery",
+                ),
+            )
 
     for t, v in self.water_balance.items():
         iscale.constraint_scaling_transform(
@@ -320,6 +321,15 @@ def calculate_scaling_factors_sido(self):
             ),
         )  # would just be a duplicate of above
 
+    if hasattr(self, "eq_isobaric"):
+        for (t, port), condata in self.eq_isobaric.items():
+            obj = getattr(self, port)
+            sf_P = min(
+                iscale.get_scaling_factor(self.inlet.pressure[t], default=1),
+                iscale.get_scaling_factor(obj.pressure[t], default=1),
+            )
+            iscale.constraint_scaling_transform(condata, sf_P)
+
 
 def _get_Q_sido(self, t):
     return self.properties_in[t].flow_vol

watertap/flowsheets/full_water_resource_recovery_facility/test/test_full_WRRF_with_ASM1_ADM1.py

@@ -230,7 +230,8 @@ def test_optimization_windows(self, optimized_system_frame):
             # Python 3.9 and 3.10
             cond == pytest.approx(1.95367e11, rel=1e-2)
             # Python 3.11 and 3.12
-            or cond == pytest.approx(3.44132e11, rel=1e-2)
+            # or cond == pytest.approx(3.44132e11, rel=1e-2)
+            or cond == pytest.approx(2.71713e11, rel=1e-2)
         )
 
     @pytest.mark.requires_idaes_solver
@@ -245,5 +246,7 @@ def test_optimization_linux(self, optimized_system_frame):
         # Check condition number to confirm scaling
         jac, _ = get_jacobian(m.rescaled_model, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            3.44152e11, rel=1e-3
+            3.44152e11,
+            # 2.71713e11,
+            rel=1e-3,
         )

watertap/property_models/multicomp_aq_sol_prop_pack.py

@@ -1390,6 +1390,28 @@ def build(self):
             doc="State pressure",
         )
 
+        if self.params.config.material_flow_basis == MaterialFlowBasis.mass:
+            self.flow_mass_phase_comp = Var(
+                self.params.phase_list,
+                self.params.component_list,
+                initialize=0.5,
+                bounds=(0, None),
+                units=pyunits.kg / pyunits.s,
+                doc="Component Mass flowrate",
+            )
+        elif self.params.config.material_flow_basis == MaterialFlowBasis.molar:
+            self.flow_mol_phase_comp = Var(
+                self.params.phase_list,
+                self.params.component_list,
+                initialize=0.1,
+                bounds=(0, None),
+                domain=NonNegativeReals,
+                units=pyunits.mol / pyunits.s,
+                doc="Component molar flow rate",
+            )
+        else:
+            raise ConfigurationError()
+
     # -----------------------------------------------------------------------------
     # Property Methods
     # Material flow state variables generated via on-demand props
@@ -1407,8 +1429,8 @@ def _flow_mol_phase_comp(self):
 
             def rule_flow_mol_phase_comp(b, p, j):
                 return (
-                    b.flow_mass_phase_comp[p, j]
-                    == b.flow_mol_phase_comp[p, j] * b.params.mw_comp[j]
+                    b.flow_mass_phase_comp[p, j] / b.params.mw_comp[j]
+                    == b.flow_mol_phase_comp[p, j]
                 )
 
             self.eq_flow_mol_phase_comp = Constraint(
@@ -1441,15 +1463,9 @@ def rule_flow_mass_phase_comp(b, p, j):
             )
 
     def _flow_mass_comp(self):
-        add_object_reference(
-            self,
-            "flow_mass_comp",
-            {
-                j: self.flow_mass_phase_comp[p, j]
-                for j in self.params.component_list
-                for p in self.params.phase_list
-            },
-        )
+        @self.Expression(self.params.component_list)
+        def flow_mass_comp(b, j):
+            return b.flow_mass_phase_comp["Liq", j]
 
     def _mass_frac_phase_comp(self):
         self.mass_frac_phase_comp = Var(
@@ -2781,6 +2797,22 @@ def calculate_scaling_factors(self):
                             self.flow_mol_phase_comp["Liq", j], sf
                         )
 
+        if self.is_property_constructed("flow_mol_comp"):
+            for j in self.flow_mol_comp:
+                iscale.set_scaling_factor(
+                    self.flow_mol_comp[j],
+                    # Don't provide a default since this should always be set
+                    iscale.get_scaling_factor(self.flow_mol_phase_comp["Liq", j]),
+                )
+
+        if self.is_property_constructed("flow_mass_comp"):
+            for j in self.flow_mass_comp:
+                iscale.set_scaling_factor(
+                    self.flow_mass_comp[j],
+                    # Don't provide a default since this should always be set
+                    iscale.get_scaling_factor(self.flow_mass_phase_comp["Liq", j]),
+                )
+
         # The following variables and parameters have computed scaling factors;
         # Users do not have to input scaling factors but, if they do, their value
         # will override.

watertap/property_models/seawater_prop_pack.py

@@ -754,6 +754,7 @@ def define_metadata(cls, obj):
                 "molality_phase_comp": {"method": "_molality_phase_comp"},
                 "visc_d_phase": {"method": "_visc_d_phase"},
                 "pressure_osm_phase": {"method": "_pressure_osm_phase"},
+                "energy_density_phase": {"method": "_energy_density_phase"},
                 "enth_mass_phase": {"method": "_enth_mass_phase"},
                 "pressure_sat": {"method": "_pressure_sat"},
                 "cp_mass_phase": {"method": "_cp_mass_phase"},
@@ -1419,6 +1420,17 @@ def rule_enth_mass_phase(b, p):
             self.params.phase_list, rule=rule_enth_mass_phase
         )
 
+    def _energy_density_phase(self):
+        def rule_energy_density_phase(b, p):
+            # We have b.energy_inernal_mass_phase
+            # = b.enth_mass_phase[p] - b.pressure/b.dens_mass_phase[p]
+            # Energy density is then just b.energy_inernal_mass_phase * b.dens_mass_phase[p]
+            return self.enth_mass_phase[p] * self.dens_mass_phase[p] - b.pressure
+
+        self.energy_density_phase = Expression(
+            self.params.phase_list, rule=rule_energy_density_phase
+        )
+
     def _enth_flow(self):
         # enthalpy flow expression for get_enthalpy_flow_terms method
 
@@ -1597,12 +1609,13 @@ def get_enthalpy_flow_terms(self, p):
         """Create enthalpy flow terms."""
         return self.enth_flow
 
-    # TODO: make property package compatible with dynamics
-    # def get_material_density_terms(self, p, j):
-    #     """Create material density terms."""
+    def get_material_density_terms(self, p, j):
+        """Create material density terms."""
+        return self.conc_mass_phase_comp[p, j]
 
-    # def get_enthalpy_density_terms(self, p):
-    #     """Create enthalpy density terms."""
+    def get_energy_density_terms(self, p):
+        """Create energy density terms."""
+        return self.energy_density_phase[p]
 
     def default_material_balance_type(self):
         return MaterialBalanceType.componentTotal