More efficient method for check_parallel_jacobian

idaes/core/base/control_volume1d.py

@@ -2554,3 +2554,27 @@ def calculate_scaling_factors(self):
                     iscale.get_scaling_factor(self.element_accumulation[t, x, e]),
                     overwrite=False,
                 )
+
+        # Collocation (Lagrange-Legendre) support
+        if hasattr(self, "_flow_terms_length_domain_cont_eq"):
+            for (t, x, p, j), c in self._flow_terms_length_domain_cont_eq.items():
+                sf = iscale.get_scaling_factor(self._flow_terms[t, x, p, j])
+                iscale.constraint_scaling_transform(c, sf, overwrite=False)
+
+        if hasattr(self, "elemental_flow_term_length_domain_cont_eq"):
+            for (t, x, e), c in self.elemental_flow_term_length_domain_cont_eq.items():
+                sf = iscale.get_scaling_factor(self.elemental_flow_term[t, x, e])
+                iscale.constraint_scaling_transform(c, sf, overwrite=False)
+
+        if hasattr(self, "pressure_length_domain_cont_eq"):
+            for (t, x), c in self.pressure_length_domain_cont_eq.items():
+                sf_P = iscale.get_scaling_factor(
+                    self.properties[t, x].pressure, default=1, warning=True
+                )
+                iscale.constraint_scaling_transform(c, sf_P, overwrite=False)
+
+        if hasattr(self, "_enthalpy_flow_length_domain_cont_eq"):
+            for (t, x, p), c in self._enthalpy_flow_length_domain_cont_eq.items():
+                # No need for default because it should already have a scaling factor
+                sf = iscale.get_scaling_factor(self._enthalpy_flow[t, x, p])
+                iscale.constraint_scaling_transform(c, sf, overwrite=False)

idaes/core/util/model_diagnostics.py

@@ -28,7 +28,7 @@
 import numpy as np
 from scipy.linalg import svd
 from scipy.sparse.linalg import svds, norm
-from scipy.sparse import issparse, find
+from scipy.sparse import issparse, find, triu, diags
 
 from pyomo.environ import (
     Binary,
@@ -289,9 +289,11 @@ def svd_sparse(jacobian, number_singular_values):
 CONFIG.declare(
     "parallel_component_tolerance",
     ConfigValue(
-        default=1e-8,
+        default=1e-15,
         domain=float,
         description="Tolerance for identifying near-parallel Jacobian rows/columns",
+        doc="Absolute tolerance for considering two Jacobian rows/columns to be considered. "
+        "parallel. A smaller value means more stringent requirements for colinearity.",
     ),
 )
 CONFIG.declare(
@@ -3712,8 +3714,9 @@ def ipopt_solve_halt_on_error(model, options=None):
 
 def check_parallel_jacobian(
     model,
-    tolerance: float = 1e-4,
+    tolerance: float = 1e-8,
     direction: str = "row",
+    zero_norm_tolerance: float = 1e-8,
     jac=None,
     nlp=None,
 ):
@@ -3746,8 +3749,7 @@ def check_parallel_jacobian(
         list of 2-tuples containing parallel Pyomo components
 
     """
-    # Thanks to Robby Parker for the sparse matrix implementation and
-    # significant performance improvements
+    # Thanks to Robby Parker and Doug Allan for significant performance improvements
 
     if direction not in ["row", "column"]:
         raise ValueError(
@@ -3762,52 +3764,64 @@ def check_parallel_jacobian(
     # they correspond to.
     if direction == "row":
         components = nlp.get_pyomo_constraints()
-        csrjac = jac.tocsr()
-        # Make everything a column vector (CSC) for consistency
-        vectors = [csrjac[i, :].transpose().tocsc() for i in range(len(components))]
+        mat = jac.tocsr()
+
     elif direction == "column":
         components = nlp.get_pyomo_variables()
-        cscjac = jac.tocsc()
-        vectors = [cscjac[:, i] for i in range(len(components))]
-
-    # List to store pairs of parallel components
-    parallel = []
-
-    vectors_by_nz = {}
-    for vecidx, vec in enumerate(vectors):
-        maxval = max(np.abs(vec.data))
-        # Construct tuple of sorted col/row indices that participate
-        # in this vector (with non-negligible coefficient).
-        nz = tuple(
-            sorted(
-                idx
-                for idx, val in zip(vec.indices, vec.data)
-                if abs(val) > tolerance and abs(val) / maxval > tolerance
-            )
-        )
-        if nz in vectors_by_nz:
-            # Store the index as well so we know what component this
-            # correrponds to.
-            vectors_by_nz[nz].append((vec, vecidx))
-        else:
-            vectors_by_nz[nz] = [(vec, vecidx)]
-
-    for vecs in vectors_by_nz.values():
-        for idx, (u, uidx) in enumerate(vecs):
-            # idx is the "local index", uidx is the "global index"
-            # Frobenius norm of the matrix is 2-norm of this column vector
-            unorm = norm(u, ord="fro")
-            for v, vidx in vecs[idx + 1 :]:
-                vnorm = norm(v, ord="fro")
-
-                # Explicitly multiply a row vector * column vector
-                prod = u.transpose().dot(v)
-                absprod = abs(prod[0, 0])
-                diff = abs(absprod - unorm * vnorm)
-                if diff <= tolerance or diff <= tolerance * max(unorm, vnorm):
-                    parallel.append((uidx, vidx))
-
-    parallel = [(components[uidx], components[vidx]) for uidx, vidx in parallel]
+        mat = jac.transpose().tocsr()
+
+    # Take product of all rows/columns with all rows/columns by taking outer
+    # product of matrix with itself
+    outer = mat @ mat.transpose()
+
+    # Along the diagonal of the outer product you get the dot product of the row
+    # with itself, which is equal to the norm squared.
+    norms = np.sqrt(outer.diagonal())
+
+    # Want to ignore indices with zero norm. By zeroing out the corresponding
+    # entries of the scaling matrix, we set the corresponding rows and columns
+    # to zero, which will then be ignored.
+
+    zero_norm_indices = np.nonzero(np.abs(norms) <= zero_norm_tolerance)
+    inv_norms = 1 / norms
+    inv_norms[zero_norm_indices] = 0
+
+    # Divide each row and each column by the vector norm. This leaves
+    # the entries as dot(u, v) / (norm(u) * norm(v)). The exception is
+    # entries with "zero norm", whose corresponding rows and columns are
+    # set to zero.
+    scaling = diags(inv_norms)
+    outer = scaling * outer * scaling
+
+    # Get rid of duplicate values by only taking upper triangular part of
+    # resulting matrix
+    upper_tri = triu(outer)
+
+    # Set diagonal elements to zero
+    # Subtracting diags(upper_tri.diagonal()) is a more reliable
+    # method of getting the entries to exactly zero than subtracting
+    # an identity matrix, where one can encounter values of 1e-16
+    upper_tri = upper_tri - diags(upper_tri.diagonal())
+
+    # Get the nonzero entries of upper_tri in three arrays,
+    # corresponding to row indices, column indices, and values
+    rows, columns, values = find(upper_tri)
+
+    # We have that dot(u,v) == norm(u) * norm(v) * cos(theta) in which
+    # theta is the angle between u and v. If theta is approximately
+    # 0 or pi, sqrt(2*(1 - abs(dot(u,v)) / (norm(u) * norm(v)))) approximately
+    # equals the number of radians from 0 or pi. A tolerance of 1e-8 corresponds
+    # to a tolerance of sqrt(2) * 1e-4 radians
+
+    # The expression (1 - abs(values) < tolerance) returns an array with values
+    # of ones and zeros, depending on whether the condition is fulfilled or not.
+    # We then find indices where it is filled using np.nonzero.
+    parallel_1D = np.nonzero(1 - abs(values) < tolerance)[0]
+
+    parallel = [
+        (components[rows[idx]], components[columns[idx]]) for idx in parallel_1D
+    ]
+
     return parallel
 
 

idaes/core/util/tests/test_model_diagnostics.py

@@ -990,18 +990,12 @@ def test_display_near_parallel_variables(self):
 
         stream = StringIO()
         dt.display_near_parallel_variables(stream)
-
-        expected = """====================================================================================
-The following pairs of variables are nearly parallel:
-
-    v1, v2
-    v1, v4
-    v2, v4
-
-====================================================================================
-"""
-
-        assert stream.getvalue() == expected
+        expected_pairs = ["v1, v2", "v1, v4", "v2, v4"]
+        assert (
+            "The following pairs of variables are nearly parallel:" in stream.getvalue()
+        )
+        for pair in expected_pairs:
+            assert pair in stream.getvalue()
 
     @pytest.mark.component
     def test_collect_structural_warnings_base_case(self, model):
@@ -1469,7 +1463,7 @@ def test_report_numerical_issues_jacobian(self):
     WARNING: 1 Constraint with large residuals (>1.0E-05)
     WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)
     WARNING: 1 Constraint with extreme Jacobian values (<1.0E-08 or >1.0E+08)
-    WARNING: 3 pairs of variables are parallel (to tolerance 1.0E-08)
+    WARNING: 1 pair of variables are parallel (to tolerance 1.0E-08)
 
 ------------------------------------------------------------------------------------
 4 Cautions