Add testing for ScalingProfiler

watertap/unit_models/tests/test_clarifier.py

@@ -13,6 +13,7 @@
 Tests for clarifier.
 """
 __author__ = "Chenyu Wang"
+from io import StringIO
 import pytest
 from pyomo.environ import (
     ConcreteModel,
@@ -28,7 +29,7 @@
     get_jacobian,
     jacobian_cond,
 )
-
+from idaes.core.scaling.scaler_profiling import ScalingProfiler
 from watertap.core.solvers import get_solver
 
 from watertap.unit_models.tests.unit_test_harness import UnitTestHarness
@@ -625,3 +626,165 @@ def test_example_case_scaler(self):
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
             2.0028333e4, rel=1e-3
         )
+
+
+def build_model():
+    m = ConcreteModel()
+    m.fs = FlowsheetBlock(dynamic=False)
+
+    m.fs.properties = ASM1ParameterBlock()
+
+    m.fs.unit = Clarifier(
+        property_package=m.fs.properties,
+        outlet_list=["underflow", "effluent"],
+        split_basis=SplittingType.componentFlow,
+    )
+
+    m.fs.unit.inlet.temperature.fix(298.15 * units.K)
+    m.fs.unit.inlet.pressure.fix(1 * units.atm)
+    m.fs.unit.inlet.flow_vol.fix(18446 * units.m**3 / units.day)
+
+    m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "S_S"].fix(58 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(92 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(363 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "X_BH"].fix(50 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "X_BA"].fix(1e-3 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "X_P"].fix(1e-3 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "S_O"].fix(1e-3 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "S_NO"].fix(1e-3 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "S_NH"].fix(23 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "S_ND"].fix(5 * units.g / units.m**3)
+    m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(16 * units.g / units.m**3)
+
+    # Alkalinity was given in mg/L based on C
+    m.fs.unit.inlet.alkalinity[0].fix(7 * units.mol / units.m**3)
+
+    # Unit option
+    m.fs.unit.split_fraction[0, "effluent", "H2O"].fix(0.993)
+    m.fs.unit.split_fraction[0, "effluent", "S_I"].fix(0.993)
+    m.fs.unit.split_fraction[0, "effluent", "S_S"].fix(0.993)
+    m.fs.unit.split_fraction[0, "effluent", "X_I"].fix(0.5192)
+    m.fs.unit.split_fraction[0, "effluent", "X_S"].fix(0.5192)
+    m.fs.unit.split_fraction[0, "effluent", "X_BH"].fix(0.5192)
+    m.fs.unit.split_fraction[0, "effluent", "X_BA"].fix(0.5192)
+    m.fs.unit.split_fraction[0, "effluent", "X_P"].fix(0.5192)
+    m.fs.unit.split_fraction[0, "effluent", "S_O"].fix(0.993)
+    m.fs.unit.split_fraction[0, "effluent", "S_NO"].fix(0.993)
+    m.fs.unit.split_fraction[0, "effluent", "S_NH"].fix(0.993)
+    m.fs.unit.split_fraction[0, "effluent", "S_ND"].fix(0.993)
+    m.fs.unit.split_fraction[0, "effluent", "X_ND"].fix(0.5192)
+    m.fs.unit.split_fraction[0, "effluent", "S_ALK"].fix(0.993)
+
+    solver = get_solver()
+    solver.solve(m)
+
+    return m
+
+
+def scale_vars_with_scalers(m):
+    scaler = ClarifierScaler()
+    scaler.scale_model(
+        m.fs.unit,
+        submodel_scalers={
+            m.fs.unit.mixed_state: ASM1PropertiesScaler,
+            m.fs.unit.underflow_state: ASM1PropertiesScaler,
+            m.fs.unit.effluent_state: ASM1PropertiesScaler,
+        },
+    )
+
+
+def scale_vars_with_iscale(m):
+    # Set scaling factors for badly scaled variables
+    iscale.set_scaling_factor(m.fs.unit.underflow_state[0.0].pressure, 1e-5)
+    iscale.set_scaling_factor(
+        m.fs.unit.underflow_state[0.0].conc_mass_comp["X_BA"], 1e3
+    )
+    iscale.set_scaling_factor(m.fs.unit.underflow_state[0.0].conc_mass_comp["X_P"], 1e3)
+    iscale.set_scaling_factor(m.fs.unit.underflow_state[0.0].conc_mass_comp["S_O"], 1e3)
+    iscale.set_scaling_factor(
+        m.fs.unit.underflow_state[0.0].conc_mass_comp["S_NO"], 1e3
+    )
+    iscale.set_scaling_factor(m.fs.unit.effluent_state[0.0].pressure, 1e-5)
+    iscale.set_scaling_factor(m.fs.unit.effluent_state[0.0].conc_mass_comp["X_BA"], 1e7)
+    iscale.set_scaling_factor(m.fs.unit.effluent_state[0.0].conc_mass_comp["X_P"], 1e7)
+    iscale.set_scaling_factor(m.fs.unit.effluent_state[0.0].conc_mass_comp["S_O"], 1e7)
+    iscale.set_scaling_factor(m.fs.unit.effluent_state[0.0].conc_mass_comp["S_NO"], 1e7)
+
+    iscale.calculate_scaling_factors(m.fs.unit)
+
+
+def perturb_solution(m):
+    m.fs.unit.inlet.flow_vol.fix(18446 * 0.8 * units.m**3 / units.day)
+    m.fs.unit.split_fraction[0, "effluent", "H2O"].fix(0.993 * 0.35)
+    m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * 1.5 * units.g / units.m**3)
+
+
+@pytest.mark.requires_idaes_solver
+@pytest.mark.unit
+def test_scaling_profiler_with_scalers():
+    sp = ScalingProfiler(
+        build_model=build_model,
+        user_scaling=scale_vars_with_scalers,
+        perturb_state=perturb_solution,
+    )
+
+    stream = StringIO()
+
+    sp.report_scaling_profiles(stream=stream)
+
+    expected = """
+============================================================================
+Scaling Profile Report
+----------------------------------------------------------------------------
+Scaling Method           || User Scaling           || Perfect Scaling
+Unscaled                 || 6.241E+06 | Solved 3   ||
+Vars Only                || 1.359E+10 | Solved 3   || 1.356E+14 | Solved 1  
+Harmonic                 || 1.359E+10 | Solved 3   || 1.228E+02 | Solved 1  
+Inverse Sum              || 1.359E+10 | Solved 3   || 6.636E+03 | Solved 1  
+Inverse Root Sum Squares || 1.359E+10 | Solved 3   || 6.633E+03 | Solved 1  
+Inverse Maximum          || 1.359E+10 | Solved 3   || 6.636E+03 | Solved 1  
+Inverse Minimum          || 1.359E+10 | Solved 3   || 9.327E+01 | Solved 1  
+Nominal L1 Norm          || 1.359E+10 | Solved 3   || 1.817E+03 | Solved 1  
+Nominal L2 Norm          || 1.359E+10 | Solved 3   || 3.354E+03 | Solved 1  
+Actual L1 Norm           || 1.359E+10 | Solved 3   || 9.450E+01 | Solved 1  
+Actual L2 Norm           || 1.359E+10 | Solved 3   || 8.480E+01 | Solved 1  
+============================================================================
+"""
+
+    assert stream.getvalue() == expected
+
+
+@pytest.mark.requires_idaes_solver
+@pytest.mark.unit
+def test_scaling_profiler_with_iscale():
+    sp = ScalingProfiler(
+        build_model=build_model,
+        user_scaling=scale_vars_with_iscale,
+        perturb_state=perturb_solution,
+    )
+
+    stream = StringIO()
+
+    sp.report_scaling_profiles(stream=stream)
+
+    expected = """
+============================================================================
+Scaling Profile Report
+----------------------------------------------------------------------------
+Scaling Method           || User Scaling           || Perfect Scaling
+Unscaled                 || 6.241E+06 | Solved 3   ||
+Vars Only                || 2.999E+09 | Solved 2   || 1.356E+14 | Solved 1  
+Harmonic                 || 1.027E+06 | Solved 2   || 1.228E+02 | Solved 1  
+Inverse Sum              || 8.148E+06 | Solved 2   || 6.636E+03 | Solved 1  
+Inverse Root Sum Squares || 8.114E+06 | Solved 2   || 6.633E+03 | Solved 1  
+Inverse Maximum          || 8.139E+06 | Solved 2   || 6.636E+03 | Solved 1  
+Inverse Minimum          || 1.257E+06 | Solved 2   || 9.327E+01 | Solved 1  
+Nominal L1 Norm          || 1.632E+09 | Solved 2   || 1.817E+03 | Solved 1  
+Nominal L2 Norm          || 1.972E+09 | Solved 2   || 3.354E+03 | Solved 1  
+Actual L1 Norm           || 1.776E+05 | Solved 2   || 9.450E+01 | Solved 1  
+Actual L2 Norm           || 1.723E+05 | Solved 2   || 8.480E+01 | Solved 1  
+============================================================================
+"""
+
+    assert stream.getvalue() == expected