Resolving Potential Evaluation Errors

watertap/flowsheets/dye_desalination/tests/test_dye_desalination.py

@@ -303,13 +303,13 @@ def test_solve(self, system_frame):
             m.fs.brine.flow_mass_phase_comp[0, "Liq", "H2O"]
         )
 
-        assert pytest.approx(-45144.537, rel=1e-5) == value(
+        assert pytest.approx(-43659.193, rel=1e-5) == value(
             m.fs.desalination.RO.deltaP[0]
         )
 
-        assert pytest.approx(1299.7286, rel=1e-5) == value(m.fs.desalination.RO.area)
+        assert pytest.approx(1256.9627, rel=1e-5) == value(m.fs.desalination.RO.area)
 
-        assert pytest.approx(3.26732341e-7, rel=1e-5) == value(
+        assert pytest.approx(1.649766198e-7, rel=1e-5) == value(
             m.fs.desalination.RO.flux_mass_phase_comp[0, 1, "Liq", "TDS"]
         )
 

watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py

@@ -556,7 +556,7 @@ def scale_variables(m):
         )
         iscale.set_scaling_factor(block.control_volume.material_balances, 1e3)
 
-    iscale.set_scaling_factor(m.fs.AD.KH_co2, 1e1)
+    iscale.set_scaling_factor(m.fs.AD.KH_co2, 1e2)
     iscale.set_scaling_factor(m.fs.AD.KH_ch4, 1e1)
     iscale.set_scaling_factor(m.fs.AD.KH_h2, 1e2)
 

watertap/property_models/NaCl_T_dep_prop_pack.py

@@ -1026,7 +1026,7 @@ def _mass_frac_phase_comp(self):
         )
 
         def rule_mass_frac_phase_comp(b, p, j):
-            return b.mass_frac_phase_comp[p, j] == b.flow_mass_phase_comp[p, j] / sum(
+            return b.flow_mass_phase_comp[p, j] == b.mass_frac_phase_comp[p, j] * sum(
                 b.flow_mass_phase_comp[p, j] for j in self.params.component_list
             )
 
@@ -1160,7 +1160,7 @@ def _mole_frac_phase_comp(self):
         )
 
         def rule_mole_frac_phase_comp(b, p, j):
-            return b.mole_frac_phase_comp[p, j] == b.flow_mol_phase_comp[p, j] / sum(
+            return b.flow_mol_phase_comp[p, j] == b.mole_frac_phase_comp[p, j] * sum(
                 b.flow_mol_phase_comp[p, j] for j in b.params.component_list
             )
 

watertap/property_models/NaCl_prop_pack.py

@@ -512,7 +512,7 @@ def _mass_frac_phase_comp(self):
         )
 
         def rule_mass_frac_phase_comp(b, p, j):
-            return b.mass_frac_phase_comp[p, j] == b.flow_mass_phase_comp[p, j] / sum(
+            return b.flow_mass_phase_comp[p, j] == b.mass_frac_phase_comp[p, j] * sum(
                 b.flow_mass_phase_comp[p, j] for j in self.params.component_list
             )
 
@@ -623,7 +623,7 @@ def _mole_frac_phase_comp(self):
         )
 
         def rule_mole_frac_phase_comp(b, p, j):
-            return b.mole_frac_phase_comp[p, j] == b.flow_mol_phase_comp[p, j] / sum(
+            return b.flow_mol_phase_comp[p, j] == b.mole_frac_phase_comp[p, j] * sum(
                 b.flow_mol_phase_comp[p, j] for j in b.params.component_list
             )
 

watertap/property_models/multicomp_aq_sol_prop_pack.py

@@ -1181,7 +1181,7 @@ def _mass_frac_phase_comp(self):
         )
 
         def rule_mass_frac_phase_comp(b, p, j):
-            return b.mass_frac_phase_comp[p, j] == b.flow_mass_phase_comp[p, j] / sum(
+            return b.flow_mass_phase_comp[p, j] == b.mass_frac_phase_comp[p, j] * sum(
                 b.flow_mass_phase_comp[p, j] for j in self.params.component_list
             )
 
@@ -1370,7 +1370,7 @@ def _mole_frac_phase_comp(self):
         )
 
         def rule_mole_frac_phase_comp(b, p, j):
-            return b.mole_frac_phase_comp[p, j] == b.flow_mol_phase_comp[p, j] / sum(
+            return b.flow_mol_phase_comp[p, j] == b.mole_frac_phase_comp[p, j] * sum(
                 b.flow_mol_phase_comp[p, j] for j in b.params.component_list
             )
 

watertap/property_models/seawater_prop_pack.py

@@ -1074,7 +1074,7 @@ def _mass_frac_phase_comp(self):
         )
 
         def rule_mass_frac_phase_comp(b, p, j):
-            return b.mass_frac_phase_comp[p, j] == b.flow_mass_phase_comp[p, j] / sum(
+            return b.flow_mass_phase_comp[p, j] == b.mass_frac_phase_comp[p, j] * sum(
                 b.flow_mass_phase_comp[p, j] for j in b.params.component_list
             )
 
@@ -1214,7 +1214,7 @@ def _mole_frac_phase_comp(self):
         )
 
         def rule_mole_frac_phase_comp(b, p, j):
-            return b.mole_frac_phase_comp[p, j] == b.flow_mol_phase_comp[p, j] / sum(
+            return b.flow_mol_phase_comp[p, j] == b.mole_frac_phase_comp[p, j] * sum(
                 b.flow_mol_phase_comp[p, j] for j in b.params.component_list
             )
 

watertap/property_models/unit_specific/NDMA_prop_pack.py

@@ -448,7 +448,7 @@ def _mass_frac_phase_comp(self):
         )
 
         def rule_mass_frac_phase_comp(b, p, j):
-            return b.mass_frac_phase_comp[p, j] == b.flow_mass_phase_comp[p, j] / sum(
+            return b.flow_mass_phase_comp[p, j] == b.mass_frac_phase_comp[p, j] * sum(
                 b.flow_mass_phase_comp[p, j] for j in self.params.component_list
             )
 
@@ -556,7 +556,7 @@ def _mole_frac_phase_comp(self):
         )
 
         def rule_mole_frac_phase_comp(b, p, j):
-            return b.mole_frac_phase_comp[p, j] == b.flow_mol_phase_comp[p, j] / sum(
+            return b.flow_mol_phase_comp[p, j] == b.mole_frac_phase_comp[p, j] * sum(
                 b.flow_mol_phase_comp[p, j] for j in b.params.component_list
             )
 

watertap/property_models/unit_specific/activated_sludge/asm1_reactions.py

@@ -413,6 +413,7 @@ def _rxn_rate(self):
         try:
 
             def rate_expression_rule(b, r):
+                eps = 1e-30 * pyo.units.kg / pyo.units.m**3
                 if r == "R1":
                     # R1: Aerobic growth of heterotrophs
                     return b.reaction_rate[r] == pyo.units.convert(
@@ -491,6 +492,7 @@ def rate_expression_rule(b, r):
                         / (
                             b.params.K_X * b.conc_mass_comp_ref["X_BH"]
                             + b.conc_mass_comp_ref["X_S"]
+                            + eps
                         )
                         * (
                             (

watertap/property_models/unit_specific/activated_sludge/modified_asm2d_reactions.py

@@ -1056,7 +1056,6 @@ def build(self):
         super().build()
 
         # Create references to state vars
-        # Concentration
         add_object_reference(self, "conc_mass_comp_ref", self.state_ref.conc_mass_comp)
 
     # Rate of reaction method
@@ -1071,6 +1070,7 @@ def _rxn_rate(self):
         try:
             # TODO: Refer to Flores-Alsina c code to account for sulfur in rate expressions
             def rate_expression_rule(b, r):
+                eps = 1e-30 * pyo.units.kg / pyo.units.m**3
                 if r == "R1":
                     # R1: Aerobic hydrolysis
                     return b.reaction_rate[r] == pyo.units.convert(
@@ -1084,6 +1084,7 @@ def rate_expression_rule(b, r):
                             / (
                                 b.params.KL_X * b.conc_mass_comp_ref["X_H"]
                                 + b.conc_mass_comp_ref["X_S"]
+                                + eps
                             )
                         )
                         * b.conc_mass_comp_ref["X_H"],
@@ -1107,6 +1108,7 @@ def rate_expression_rule(b, r):
                             / (
                                 b.params.KL_X * b.conc_mass_comp_ref["X_H"]
                                 + b.conc_mass_comp_ref["X_S"]
+                                + eps
                             )
                         )
                         * b.conc_mass_comp_ref["X_H"],
@@ -1130,6 +1132,7 @@ def rate_expression_rule(b, r):
                             / (
                                 b.params.KL_X * b.conc_mass_comp_ref["X_H"]
                                 + b.conc_mass_comp_ref["X_S"]
+                                + eps
                             )
                         )
                         * b.conc_mass_comp_ref["X_H"],
@@ -1323,6 +1326,7 @@ def rate_expression_rule(b, r):
                             / (
                                 b.params.KP_PP * b.conc_mass_comp_ref["X_PAO"]
                                 + b.conc_mass_comp_ref["X_PP"]
+                                + eps
                             )
                         )
                         * b.conc_mass_comp_ref["X_PAO"],
@@ -1345,6 +1349,7 @@ def rate_expression_rule(b, r):
                             / (
                                 b.params.KP_PHA * b.conc_mass_comp_ref["X_PAO"]
                                 + b.conc_mass_comp_ref["X_PHA"]
+                                + eps
                             )
                         )
                         * (
@@ -1383,6 +1388,7 @@ def rate_expression_rule(b, r):
                             / (
                                 b.params.KP_PHA * b.conc_mass_comp_ref["X_PAO"]
                                 + b.conc_mass_comp_ref["X_PHA"]
+                                + eps
                             )
                         )
                         * (
@@ -1394,6 +1400,7 @@ def rate_expression_rule(b, r):
                                 b.params.KI_PP * b.conc_mass_comp_ref["X_PAO"]
                                 + b.params.K_MAX * b.conc_mass_comp_ref["X_PAO"]
                                 - b.conc_mass_comp_ref["X_PP"]
+                                + eps
                             )
                         )
                         * b.conc_mass_comp_ref["X_PAO"],
@@ -1420,6 +1427,7 @@ def rate_expression_rule(b, r):
                             / (
                                 b.params.KP_PHA * b.conc_mass_comp_ref["X_PAO"]
                                 + b.conc_mass_comp_ref["X_PHA"]
+                                + eps
                             )
                         )
                         * b.conc_mass_comp_ref["X_PAO"],
@@ -1451,6 +1459,7 @@ def rate_expression_rule(b, r):
                             / (
                                 b.params.KP_PHA * b.conc_mass_comp_ref["X_PAO"]
                                 + b.conc_mass_comp_ref["X_PHA"]
+                                + eps
                             )
                         )
                         * b.conc_mass_comp_ref["X_PAO"],

watertap/property_models/unit_specific/anaerobic_digestion/adm1_properties_vapor.py

@@ -272,6 +272,8 @@ def build(self):
             units=pyo.units.Pa,
         )
 
+        eps = 1e-30
+
         def pressure_sat_rule(b, j):
             if j == "S_h2":
                 return b.pressure_sat[j] == pyo.units.convert(
@@ -292,7 +294,7 @@ def pressure_sat_rule(b, j):
                     to_units=pyo.units.Pa,
                 )
             elif j == "H2O":
-                return pyo.log(b.pressure_sat[j] / pyo.units.Pa) == (
+                return pyo.log(b.pressure_sat[j] / pyo.units.Pa + eps) == (
                     pyo.log(0.0313)
                     + 5290
                     * pyo.units.K