From 62a6346a5a184b852c5d7ec3c4e1bb817b755719 Mon Sep 17 00:00:00 2001 From: Doug A Date: Wed, 10 Apr 2024 16:51:51 -0400 Subject: [PATCH 01/28] implement new method --- idaes/core/util/model_diagnostics.py | 67 +++++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 67b7adc0f4..0a44553009 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -27,7 +27,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 from pyomo.environ import ( Binary, @@ -3653,6 +3653,71 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro jac, nlp = get_jacobian(model, scaled=False) + # Get vectors that we will check, and the Pyomo components + # they correspond to. + if direction == "row": + components = nlp.get_pyomo_constraints() + mat = jac.tocsr() + + elif direction == "column": + components = nlp.get_pyomo_variables() + mat = jac.transpose().tocsr() + + norms = [norm(mat[i, :], ord="fro") for i in range(len(components))] + + # Take product of all rows/columns with all rows/columns by taking outer + # product of matrix with itself + outer = mat @ mat.transpose() + # Get rid of duplicate values by only taking upper triangular part of + # resulting matrix + upper_tri = triu(outer) + # List to store pairs of parallel components + parallel = [] + + for row, col, val in zip(*find(upper_tri), strict=True): + if row == col: + # A vector is parallel to itself + continue + diff = abs(abs(val) - norms[row] * norms[col]) + if diff <= tolerance or diff <= tolerance * max(norms[row], norms[col]): + parallel.append((components[row], components[col])) + + return parallel + +def check_parallel_jacobian_old(model, tolerance: float = 1e-4, direction: str = "row"): + """ + Check for near-parallel rows or columns in the Jacobian. + + Near-parallel rows or columns indicate a potential degeneracy in the model, + as this means that the associated constraints or variables are (near) + duplicates of each other. + + This method is based on work published in: + + Klotz, E., Identification, Assessment, and Correction of Ill-Conditioning and + Numerical Instability in Linear and Integer Programs, Informs 2014, pgs. 54-108 + https://pubsonline.informs.org/doi/epdf/10.1287/educ.2014.0130 + + Args: + model: model to be analysed + tolerance: tolerance to use to determine if constraints/variables are parallel + direction: 'row' (default, constraints) or 'column' (variables) + + Returns: + list of 2-tuples containing parallel Pyomo components + + """ + # Thanks to Robby Parker for the sparse matrix implementation and + # significant performance improvements + + if direction not in ["row", "column"]: + raise ValueError( + f"Unrecognised value for direction ({direction}). " + "Must be 'row' or 'column'." + ) + + jac, nlp = get_jacobian(model, scaled=False) + # Get vectors that we will check, and the Pyomo components # they correspond to. if direction == "row": From 6a462297b218c5bfba315eb162aaadc7727f5e33 Mon Sep 17 00:00:00 2001 From: Doug A Date: Wed, 10 Apr 2024 17:28:50 -0400 Subject: [PATCH 02/28] run black --- idaes/core/util/model_diagnostics.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 0a44553009..4f2d5f4481 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3681,9 +3681,10 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro diff = abs(abs(val) - norms[row] * norms[col]) if diff <= tolerance or diff <= tolerance * max(norms[row], norms[col]): parallel.append((components[row], components[col])) - + return parallel + def check_parallel_jacobian_old(model, tolerance: float = 1e-4, direction: str = "row"): """ Check for near-parallel rows or columns in the Jacobian. From 5e4da397ddb81c77aaa493a9ff4206cfa809949d Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 11 Apr 2024 09:11:51 -0400 Subject: [PATCH 03/28] Nonstrict zip --- idaes/core/util/model_diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 4f2d5f4481..a11bd8390c 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3674,7 +3674,7 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro # List to store pairs of parallel components parallel = [] - for row, col, val in zip(*find(upper_tri), strict=True): + for row, col, val in zip(*find(upper_tri)): if row == col: # A vector is parallel to itself continue From 8fbb339886ed0f0b94667a1941066b6dcc118d90 Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 11 Apr 2024 12:15:05 -0400 Subject: [PATCH 04/28] only relative criterion, filter out zero norm vectors --- idaes/core/util/model_diagnostics.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index a11bd8390c..fedcc232b3 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3619,7 +3619,7 @@ def ipopt_solve_halt_on_error(model, options=None): ) -def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "row"): +def check_parallel_jacobian(model, tolerance: float = 1e-8, direction: str = "row", zero_norm_tolerance: float = 1e-8): """ Check for near-parallel rows or columns in the Jacobian. @@ -3663,7 +3663,12 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro components = nlp.get_pyomo_variables() mat = jac.transpose().tocsr() - norms = [norm(mat[i, :], ord="fro") for i in range(len(components))] + norms = np.NaN(len(components)) + + for i in range(len(components)): + norms[i] = norm(mat[i, :], ord="fro") + #norms = np.array([norm(mat[i, :], ord="fro") for i in range(len(components))]) + zero_norm_indices = np.nonzero(np.abs(norms) <= zero_norm_tolerance) # Take product of all rows/columns with all rows/columns by taking outer # product of matrix with itself @@ -3678,8 +3683,16 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro if row == col: # A vector is parallel to itself continue - diff = abs(abs(val) - norms[row] * norms[col]) - if diff <= tolerance or diff <= tolerance * max(norms[row], norms[col]): + elif row in zero_norm_indices or col in zero_norm_indices: + # Any index with effectively zero norm will be reported + # by the jacobian_rows and jacobian_columns functions + continue + # 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 + elif 1 - abs(val) / (norms[row] * norms[col]) <= tolerance: parallel.append((components[row], components[col])) return parallel From da664b0a63e1fa1655ec1d5159656f50247bc16a Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 11 Apr 2024 12:20:09 -0400 Subject: [PATCH 05/28] Correct preallocation --- idaes/core/util/model_diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index fedcc232b3..9b02078baf 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3663,7 +3663,7 @@ def check_parallel_jacobian(model, tolerance: float = 1e-8, direction: str = "ro components = nlp.get_pyomo_variables() mat = jac.transpose().tocsr() - norms = np.NaN(len(components)) + norms = np.NaN * np.ones(len(components)) for i in range(len(components)): norms[i] = norm(mat[i, :], ord="fro") From 9b8a25d23362b24f2d0108ded0aa565ae5ae8957 Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 11 Apr 2024 12:31:35 -0400 Subject: [PATCH 06/28] run Black --- idaes/core/util/model_diagnostics.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 9b02078baf..32d5f136dc 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3619,7 +3619,12 @@ def ipopt_solve_halt_on_error(model, options=None): ) -def check_parallel_jacobian(model, tolerance: float = 1e-8, direction: str = "row", zero_norm_tolerance: float = 1e-8): +def check_parallel_jacobian( + model, + tolerance: float = 1e-8, + direction: str = "row", + zero_norm_tolerance: float = 1e-8, +): """ Check for near-parallel rows or columns in the Jacobian. @@ -3667,7 +3672,7 @@ def check_parallel_jacobian(model, tolerance: float = 1e-8, direction: str = "ro for i in range(len(components)): norms[i] = norm(mat[i, :], ord="fro") - #norms = np.array([norm(mat[i, :], ord="fro") for i in range(len(components))]) + # norms = np.array([norm(mat[i, :], ord="fro") for i in range(len(components))]) zero_norm_indices = np.nonzero(np.abs(norms) <= zero_norm_tolerance) # Take product of all rows/columns with all rows/columns by taking outer From 7fe9c5efc044ccfe6a40d2d30081ce368e6126d3 Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 11 Apr 2024 13:49:48 -0400 Subject: [PATCH 07/28] efficiency --- idaes/core/util/model_diagnostics.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 32d5f136dc..710df14b3b 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3668,22 +3668,27 @@ def check_parallel_jacobian( components = nlp.get_pyomo_variables() mat = jac.transpose().tocsr() - norms = np.NaN * np.ones(len(components)) + # 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 matrix + # with itself, which is equal to the norm squared. + norms = np.sqrt(outer.diagonal()) - for i in range(len(components)): - norms[i] = norm(mat[i, :], ord="fro") # norms = np.array([norm(mat[i, :], ord="fro") for i in range(len(components))]) zero_norm_indices = np.nonzero(np.abs(norms) <= zero_norm_tolerance) - # Take product of all rows/columns with all rows/columns by taking outer - # product of matrix with itself - outer = mat @ mat.transpose() # Get rid of duplicate values by only taking upper triangular part of # resulting matrix upper_tri = triu(outer) - # List to store pairs of parallel components + parallel = [] + # find returns a tuple of three arrays, with entries corresponding to matrix + # row number, matrix column number, and entry value for nonzero entries of + # upper_tri. * expands the tuple into three arguments for zip, which then + # allows us to iterate over row, column, and value simultaneously. for row, col, val in zip(*find(upper_tri)): if row == col: # A vector is parallel to itself From ab300492e9ef73cfc5e3f40c94ad7dcbb5556e07 Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 11 Apr 2024 17:27:42 -0400 Subject: [PATCH 08/28] excessive optimization --- idaes/core/util/model_diagnostics.py | 49 ++++++++++++++-------------- 1 file changed, 25 insertions(+), 24 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 710df14b3b..af2111cacb 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -27,7 +27,7 @@ import numpy as np from scipy.linalg import svd from scipy.sparse.linalg import svds, norm -from scipy.sparse import issparse, find, triu +from scipy.sparse import issparse, find, triu, diags from pyomo.environ import ( Binary, @@ -3656,7 +3656,8 @@ def check_parallel_jacobian( "Must be 'row' or 'column'." ) - jac, nlp = get_jacobian(model, scaled=False) + if nlp is None or jac is None: + jac, nlp = get_jacobian(model, scaled=False) # Get vectors that we will check, and the Pyomo components # they correspond to. @@ -3672,38 +3673,38 @@ def check_parallel_jacobian( # product of matrix with itself outer = mat @ mat.transpose() - # Along the diagonal of the outer product you get the dot product of the matrix + # 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()) - # norms = np.array([norm(mat[i, :], ord="fro") for i in range(len(components))]) zero_norm_indices = np.nonzero(np.abs(norms) <= zero_norm_tolerance) + # Want to ignore indices with zero norm. If we give them + # norms of infinity, we scale the rows and columns to zero + # which allows us to ignore them + norms[zero_norm_indices] = float('inf') # 1/float('inf') == 0 + + scaling = diags(1/norms) + outer = scaling * outer * scaling + # Get rid of duplicate values by only taking upper triangular part of # resulting matrix upper_tri = triu(outer) - parallel = [] + # 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) + + # Find values with an absolute value less than tolerance away from 1 + parallel_1D = np.nonzero(1 - abs(values) < tolerance)[0] - # find returns a tuple of three arrays, with entries corresponding to matrix - # row number, matrix column number, and entry value for nonzero entries of - # upper_tri. * expands the tuple into three arguments for zip, which then - # allows us to iterate over row, column, and value simultaneously. - for row, col, val in zip(*find(upper_tri)): - if row == col: - # A vector is parallel to itself - continue - elif row in zero_norm_indices or col in zero_norm_indices: - # Any index with effectively zero norm will be reported - # by the jacobian_rows and jacobian_columns functions - continue - # 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 - elif 1 - abs(val) / (norms[row] * norms[col]) <= tolerance: - parallel.append((components[row], components[col])) + parallel = [(components[rows[idx]], components[columns[idx]]) for idx in parallel_1D] return parallel From 8c9959c10448bad4458ed2bcf0bd5a023397dd43 Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 11 Apr 2024 17:37:55 -0400 Subject: [PATCH 09/28] add optional jacobian and nlp arguments plus fix test --- idaes/core/util/model_diagnostics.py | 2 ++ idaes/core/util/tests/test_model_diagnostics.py | 3 +++ 2 files changed, 5 insertions(+) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index af2111cacb..eb88acd767 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3624,6 +3624,8 @@ def check_parallel_jacobian( tolerance: float = 1e-8, direction: str = "row", zero_norm_tolerance: float = 1e-8, + jac = None, + nlp = None, ): """ Check for near-parallel rows or columns in the Jacobian. diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index e60c30e8b7..4ec64cc702 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -991,7 +991,10 @@ def test_display_near_parallel_variables(self): v1, v2 v1, v4 + v1, v3 v2, v4 + v2, v3 + v4, v3 ==================================================================================== """ From e4ce7071586e7d2ab8b4dc9f31f280a9be1d30d2 Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 11 Apr 2024 17:39:11 -0400 Subject: [PATCH 10/28] run black --- idaes/core/util/model_diagnostics.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index eb88acd767..600f008f2e 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3624,8 +3624,8 @@ def check_parallel_jacobian( tolerance: float = 1e-8, direction: str = "row", zero_norm_tolerance: float = 1e-8, - jac = None, - nlp = None, + jac=None, + nlp=None, ): """ Check for near-parallel rows or columns in the Jacobian. @@ -3684,9 +3684,9 @@ def check_parallel_jacobian( # Want to ignore indices with zero norm. If we give them # norms of infinity, we scale the rows and columns to zero # which allows us to ignore them - norms[zero_norm_indices] = float('inf') # 1/float('inf') == 0 + norms[zero_norm_indices] = float("inf") # 1/float('inf') == 0 - scaling = diags(1/norms) + scaling = diags(1 / norms) outer = scaling * outer * scaling # Get rid of duplicate values by only taking upper triangular part of @@ -3703,10 +3703,12 @@ def check_parallel_jacobian( # corresponding to row indices, column indices, and values rows, columns, values = find(upper_tri) - # Find values with an absolute value less than tolerance away from 1 + # Find values with an absolute value less than tolerance away from 1 parallel_1D = np.nonzero(1 - abs(values) < tolerance)[0] - parallel = [(components[rows[idx]], components[columns[idx]]) for idx in parallel_1D] + parallel = [ + (components[rows[idx]], components[columns[idx]]) for idx in parallel_1D + ] return parallel From 8438d09810ade4b60ae14a9b135660cfe382c9fb Mon Sep 17 00:00:00 2001 From: Doug A Date: Fri, 12 Apr 2024 09:22:49 -0400 Subject: [PATCH 11/28] fix failing test and add additional comments --- idaes/core/util/model_diagnostics.py | 10 ++++++++- .../core/util/tests/test_model_diagnostics.py | 21 ++++++------------- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 600f008f2e..facca4d8c4 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3686,6 +3686,10 @@ def check_parallel_jacobian( # which allows us to ignore them norms[zero_norm_indices] = float("inf") # 1/float('inf') == 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(1 / norms) outer = scaling * outer * scaling @@ -3703,7 +3707,11 @@ def check_parallel_jacobian( # corresponding to row indices, column indices, and values rows, columns, values = find(upper_tri) - # Find values with an absolute value less than tolerance away from 1 + # 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 parallel_1D = np.nonzero(1 - abs(values) < tolerance)[0] parallel = [ diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 4ec64cc702..7038e88e67 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -985,21 +985,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 - v1, v3 - v2, v4 - v2, v3 - v4, v3 - -==================================================================================== -""" - - assert stream.getvalue() == expected + expected_pairs = ["v1, v2", "v1, v4", "v1, v3", "v2, v4", "v2, v3", "v4, v3"] + 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): From e048a22bceca9a17e2ea8c88cfed31e5f786bc4c Mon Sep 17 00:00:00 2001 From: Doug A Date: Fri, 19 Apr 2024 10:30:37 -0400 Subject: [PATCH 12/28] Get rid of old method, try to make code more readable --- idaes/core/util/model_diagnostics.py | 106 +++------------------------ 1 file changed, 11 insertions(+), 95 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index facca4d8c4..344e237059 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3649,8 +3649,6 @@ 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 if direction not in ["row", "column"]: raise ValueError( @@ -3679,18 +3677,19 @@ def check_parallel_jacobian( # with itself, which is equal to the norm squared. norms = np.sqrt(outer.diagonal()) - zero_norm_indices = np.nonzero(np.abs(norms) <= zero_norm_tolerance) + # 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. - # Want to ignore indices with zero norm. If we give them - # norms of infinity, we scale the rows and columns to zero - # which allows us to ignore them - norms[zero_norm_indices] = float("inf") # 1/float('inf') == 0 + 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(1 / norms) + scaling = diags(inv_norms) outer = scaling * outer * scaling # Get rid of duplicate values by only taking upper triangular part of @@ -3712,6 +3711,10 @@ def check_parallel_jacobian( # 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 = [ @@ -3721,93 +3724,6 @@ def check_parallel_jacobian( return parallel -def check_parallel_jacobian_old(model, tolerance: float = 1e-4, direction: str = "row"): - """ - Check for near-parallel rows or columns in the Jacobian. - - Near-parallel rows or columns indicate a potential degeneracy in the model, - as this means that the associated constraints or variables are (near) - duplicates of each other. - - This method is based on work published in: - - Klotz, E., Identification, Assessment, and Correction of Ill-Conditioning and - Numerical Instability in Linear and Integer Programs, Informs 2014, pgs. 54-108 - https://pubsonline.informs.org/doi/epdf/10.1287/educ.2014.0130 - - Args: - model: model to be analysed - tolerance: tolerance to use to determine if constraints/variables are parallel - direction: 'row' (default, constraints) or 'column' (variables) - - Returns: - list of 2-tuples containing parallel Pyomo components - - """ - # Thanks to Robby Parker for the sparse matrix implementation and - # significant performance improvements - - if direction not in ["row", "column"]: - raise ValueError( - f"Unrecognised value for direction ({direction}). " - "Must be 'row' or 'column'." - ) - - jac, nlp = get_jacobian(model, scaled=False) - - # Get vectors that we will check, and the Pyomo components - # 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))] - 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] - return parallel - - def compute_ill_conditioning_certificate( model, target_feasibility_tol: float = 1e-06, From 8cedc01ba2b3e8d6a70101c29a17313b841751d9 Mon Sep 17 00:00:00 2001 From: Doug A Date: Fri, 3 May 2024 17:30:27 -0400 Subject: [PATCH 13/28] failures --- idaes/core/util/tests/test_model_diagnostics.py | 4 ++-- idaes/models/unit_models/tests/test_heat_exchanger.py | 8 +++++++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index ed4d2d6d06..16dd5512a8 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -985,7 +985,7 @@ def test_display_near_parallel_variables(self): stream = StringIO() dt.display_near_parallel_variables(stream) - expected_pairs = ["v1, v2", "v1, v4", "v1, v3", "v2, v4", "v2, v3", "v4, v3"] + expected_pairs = ["v1, v2", "v1, v4", "v2, v4"] assert ( "The following pairs of variables are nearly parallel:" in stream.getvalue() ) @@ -1450,7 +1450,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 diff --git a/idaes/models/unit_models/tests/test_heat_exchanger.py b/idaes/models/unit_models/tests/test_heat_exchanger.py index 1bcf4efd9e..85f1797fb2 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger.py @@ -23,6 +23,7 @@ ConcreteModel, Constraint, Expression, + TransformationFactory, value, Var, units as pyunits, @@ -71,6 +72,7 @@ InitializationStatus, ) from idaes.core.util import DiagnosticsToolbox +from idaes.core.util import scaling as iscale # Imports to assemble BT-PR with different units @@ -1243,7 +1245,11 @@ def test_conservation(self, iapws): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + m_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False + ) + dt = DiagnosticsToolbox(m_scaled) dt.assert_no_numerical_warnings() From 345a81397fad057b81924840b5662746b6821207 Mon Sep 17 00:00:00 2001 From: Doug A Date: Sat, 15 Jun 2024 17:08:46 -0400 Subject: [PATCH 14/28] Begun changes to address issues --- .../activity_coeff_prop_pack.py | 28 +++++--- idaes/models/unit_models/feed_flash.py | 31 ++++++++- idaes/models/unit_models/separator.py | 63 ++++++++++++++---- .../solid_liquid/tests/test_sl_separator.py | 4 +- idaes/models/unit_models/tests/test_cstr.py | 4 +- .../unit_models/tests/test_feed_flash.py | 21 +++++- idaes/models/unit_models/tests/test_flash.py | 22 ++++++- .../unit_models/tests/test_heat_exchanger.py | 65 +++++++++++++++++-- .../tests/test_heat_exchanger_1D.py | 37 +++++++++-- 9 files changed, 235 insertions(+), 40 deletions(-) diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index 6e8aab73ce..28a6eea414 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -222,7 +222,9 @@ def build(self): self.set_default_scaling("temperature_dew", 1e-1) self.set_default_scaling("temperature_bubble", 1e-1) self.set_default_scaling("flow_mol_phase", 1e-2) - self.set_default_scaling("density_mol", 1) + self.set_default_scaling("density_mol", 1/11.1e3, index="Liq") + self.set_default_scaling("density_mol", 0.31, index="Vap") + self.set_default_scaling("pressure", 1e-6) self.set_default_scaling("pressure_sat", 1e-6) self.set_default_scaling("pressure_dew", 1e-6) @@ -234,8 +236,8 @@ def build(self): self.set_default_scaling("enth_mol_phase", 1e-4, index="Vap") self.set_default_scaling("entr_mol_phase_comp", 1e-2) self.set_default_scaling("entr_mol_phase", 1e-2) - self.set_default_scaling("mw", 100) - self.set_default_scaling("mw_phase", 100) + self.set_default_scaling("mw", 1/100) + self.set_default_scaling("mw_phase", 1/100) self.set_default_scaling("gibbs_mol_phase_comp", 1e-3) self.set_default_scaling("fug_vap", 1e-6) self.set_default_scaling("fug_liq", 1e-6) @@ -2124,8 +2126,9 @@ def calculate_scaling_factors(self): if self.is_property_constructed("temperature_bubble"): iscale.set_scaling_factor(self.temperature_bubble, sf_T) if self.is_property_constructed("eq_temperature_bubble"): + sf_p = iscale.get_scaling_factor(self.pressure) iscale.constraint_scaling_transform( - self.eq_temperature_bubble, sf_T, overwrite=False + self.eq_temperature_bubble, sf_p, overwrite=False ) if self.is_property_constructed("temperature_dew"): @@ -2215,12 +2218,19 @@ def calculate_scaling_factors(self): ) iscale.constraint_scaling_transform(c, sf, overwrite=False) + if self.is_property_constructed("pressure_sat_comp"): + sf_p = iscale.get_scaling_factor(self.pressure) + for j, c in self.pressure_sat_comp.items(): + iscale.set_scaling_factor(c, sf_p, overwrite=False) + if self.is_property_constructed("eq_phase_equilibrium"): - for p, c in self.eq_phase_equilibrium.items(): - sf = iscale.get_scaling_factor( - self.eq_phase_equilibrium[p], default=1, warning=True - ) - iscale.constraint_scaling_transform(c, sf, overwrite=False) + sf_p = iscale.get_scaling_factor(self.pressure) + sf_x = { + j: iscale.get_scaling_factor(self.mole_frac_comp[j]) + for j in self.mole_frac_comp + } + for j, c in self.eq_phase_equilibrium.items(): + iscale.constraint_scaling_transform(c, sf_p*sf_x[j], overwrite=False) if self.is_property_constructed("eq_P_vap"): for p, c in self.eq_P_vap.items(): diff --git a/idaes/models/unit_models/feed_flash.py b/idaes/models/unit_models/feed_flash.py index 720d11c2b8..1cec7fefc7 100644 --- a/idaes/models/unit_models/feed_flash.py +++ b/idaes/models/unit_models/feed_flash.py @@ -16,7 +16,7 @@ from enum import Enum # Import Pyomo libraries -from pyomo.environ import Reference +from pyomo.environ import exp, log, Reference from pyomo.common.config import ConfigBlock, ConfigValue, In # Import IDAES cores @@ -31,6 +31,7 @@ from idaes.core.util.config import is_physical_parameter_block from idaes.core.util.tables import create_stream_table_dataframe from idaes.core.util.initialization import fix_state_vars +import idaes.core.util.scaling as iscale __author__ = "Andrew Lee" @@ -214,3 +215,31 @@ def fix_initialization_states(self): None """ fix_state_vars(self.control_volume.properties_in) + + def calculate_scaling_factors(self): + super().calculate_scaling_factors() + if self.config.flash_type == FlashType.isothermal: + for t in self.flowsheet().time: + sT = iscale.get_scaling_factor( + self.control_volume.properties_in[t].temperature, + default=1, + warning=True + ) + iscale.constraint_scaling_transform(self.isothermal[t], sT, overwrite=False) + elif self.config.flash_type == FlashType.isenthalpic: + cv = self.control_volume + for t in self.flowsheet().time: + s_enth = float("inf") + for p in cv.properties_in[t].phase_list: + s_enth = min( + s_enth, + iscale.get_scaling_factor( + cv.properties_in[t].get_enthalpy_flow_terms(p), + default = 1, + warning=True + ) + ) + iscale.constraint_scaling_transform(self.isenthalpic[t], s_enth, overwrite=False) + + + diff --git a/idaes/models/unit_models/separator.py b/idaes/models/unit_models/separator.py index a2a273957a..168f30aded 100644 --- a/idaes/models/unit_models/separator.py +++ b/idaes/models/unit_models/separator.py @@ -1789,26 +1789,19 @@ def calculate_scaling_factors(self): mb_type = mixed_state[t_ref].default_material_balance_type() super().calculate_scaling_factors() - if hasattr(self, "temperature_equality_eqn"): - for (t, i), c in self.temperature_equality_eqn.items(): - s = iscale.get_scaling_factor( - mixed_state[t].temperature, default=1, warning=True - ) - iscale.constraint_scaling_transform(c, s) - if hasattr(self, "pressure_equality_eqn"): for (t, i), c in self.pressure_equality_eqn.items(): s = iscale.get_scaling_factor( mixed_state[t].pressure, default=1, warning=True ) - iscale.constraint_scaling_transform(c, s) + iscale.constraint_scaling_transform(c, s, overwrite=False) if hasattr(self, "material_splitting_eqn"): if mb_type == MaterialBalanceType.componentPhase: for (t, _, p, j), c in self.material_splitting_eqn.items(): flow_term = mixed_state[t].get_material_flow_terms(p, j) s = iscale.get_scaling_factor(flow_term, default=1) - iscale.constraint_scaling_transform(c, s) + iscale.constraint_scaling_transform(c, s, overwrite=False) elif mb_type == MaterialBalanceType.componentTotal: for (t, _, j), c in self.material_splitting_eqn.items(): s = None @@ -1823,7 +1816,7 @@ def calculate_scaling_factors(self): else: _s = iscale.get_scaling_factor(ft, default=1) s = _s if _s < s else s - iscale.constraint_scaling_transform(c, s) + iscale.constraint_scaling_transform(c, s, overwrite=False) elif mb_type == MaterialBalanceType.total: pc_set = mixed_state.phase_component_set for (t, _), c in self.material_splitting_eqn.items(): @@ -1834,7 +1827,55 @@ def calculate_scaling_factors(self): else: _s = iscale.get_scaling_factor(ft, default=1) s = _s if _s < s else s - iscale.constraint_scaling_transform(c, s) + iscale.constraint_scaling_transform(c, s, overwrite=False) + + if self.config.energy_split_basis == EnergySplittingType.none: + pass + elif self.config.energy_split_basis == EnergySplittingType.equal_temperature: + for (t, i), c in self.temperature_equality_eqn.items(): + s = iscale.get_scaling_factor( + mixed_state[t].temperature, default=1, warning=True + ) + iscale.constraint_scaling_transform(c, s, overwrite=False) + + elif self.config.energy_split_basis == EnergySplittingType.equal_molar_enthalpy: + for (t, i), c in self.temperature_equality_eqn.items(): + s_enth_mol = iscale.get_scaling_factor( + mixed_state[t].enth_mol, default=1, warning=True + ) + iscale.constraint_scaling_transform(c, s_enth_mol, overwrite=False) + + elif self.config.energy_split_basis == EnergySplittingType.enthalpy_split: + # Validate split fraction type + if ( + self.config.split_basis == SplittingType.phaseComponentFlow + or self.config.split_basis == SplittingType.componentFlow + ): + raise ConfigurationError( + "{} Cannot use energy_split_basis == enthalpy_split " + "with split_basis == component or phaseComponent.".format(self.name) + ) + + for (t, i), c in self.molar_enthalpy_splitting_eqn.items(): + sf_enth = float('inf') + for p in mixed_state[t].phase_list: + sf_enth = min( + sf_enth, + iscale.get_scaling_factor( + mixed_state[t].get_enthalpy_flow_terms(p), + default=1, + warning=True + ) + ) + iscale.constraint_scaling_transform(c, sf_enth, overwrite=False) + + else: + raise BurntToast( + "{} received unrecognised value for " + "energy_split_basis. This should never happen, so" + " please contact the IDAES developers with this " + "bug.".format(self.name) + ) def _get_performance_contents(self, time_point=0): if hasattr(self, "split_fraction"): diff --git a/idaes/models/unit_models/solid_liquid/tests/test_sl_separator.py b/idaes/models/unit_models/solid_liquid/tests/test_sl_separator.py index 5ff3a5c95e..435cd3529c 100644 --- a/idaes/models/unit_models/solid_liquid/tests/test_sl_separator.py +++ b/idaes/models/unit_models/solid_liquid/tests/test_sl_separator.py @@ -194,7 +194,9 @@ def test_solve(self, model): @pytest.mark.component def test_numerical_issues(self, model): dt = DiagnosticsToolbox(model) - dt.assert_no_numerical_warnings() + # Need to scale model, perform Pyomo scaling transform, then pass that + # No scaling factors at present for Saponification + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") diff --git a/idaes/models/unit_models/tests/test_cstr.py b/idaes/models/unit_models/tests/test_cstr.py index 0701c68e32..d80e87a67c 100644 --- a/idaes/models/unit_models/tests/test_cstr.py +++ b/idaes/models/unit_models/tests/test_cstr.py @@ -247,7 +247,9 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + # Need to scale model, perform Pyomo scaling transform, then pass that + # No scaling factors at present for Saponification + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.ui @pytest.mark.unit diff --git a/idaes/models/unit_models/tests/test_feed_flash.py b/idaes/models/unit_models/tests/test_feed_flash.py index f521b2002c..23aa2b4fd1 100644 --- a/idaes/models/unit_models/tests/test_feed_flash.py +++ b/idaes/models/unit_models/tests/test_feed_flash.py @@ -20,6 +20,7 @@ from pyomo.environ import ( check_optimal_termination, ConcreteModel, + TransformationFactory, value, units as pyunits, ) @@ -44,6 +45,7 @@ InitializationStatus, ) from idaes.core.util import DiagnosticsToolbox +import idaes.core.util.scaling as iscale # ----------------------------------------------------------------------------- @@ -83,6 +85,7 @@ def btx(self): m.fs.properties = BTXParameterBlock( valid_phase=("Liq", "Vap"), activity_coeff_model="Ideal" ) + params = m.fs.properties m.fs.unit = FeedFlash(property_package=m.fs.properties) @@ -113,6 +116,10 @@ def btx(self): m.fs.unit.control_volume.properties_out[0.0].pressure_sat_comp.setlb(1e4) m.fs.unit.control_volume.properties_out[0.0].pressure_sat_comp.setub(5e6) + params.set_default_scaling("flow_mol", 1) + params.set_default_scaling("flow_mol_phase", 1) + params.set_default_scaling("flow_mol_phase_comp", 1) + return m @pytest.mark.build @@ -208,7 +215,12 @@ def test_solution(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -328,6 +340,13 @@ def test_solution(self, iapws): @pytest.mark.component def test_numerical_issues(self, iapws): dt = DiagnosticsToolbox(iapws) + + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory('core.scale_model').create_using( + iapws, + rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_flash.py b/idaes/models/unit_models/tests/test_flash.py index 7e4b459939..02e621635f 100644 --- a/idaes/models/unit_models/tests/test_flash.py +++ b/idaes/models/unit_models/tests/test_flash.py @@ -19,6 +19,7 @@ from pyomo.environ import ( check_optimal_termination, ConcreteModel, + TransformationFactory, value, units, units as pyunits, @@ -101,7 +102,7 @@ def btx(self): m.fs.properties = BTXParameterBlock( valid_phase=("Liq", "Vap"), activity_coeff_model="Ideal" ) - + m.fs.unit = Flash(property_package=m.fs.properties) m.fs.unit.inlet.flow_mol.fix(1) @@ -134,6 +135,11 @@ def btx(self): m.fs.unit.control_volume.properties_out[0.0].pressure_sat_comp.setlb(1e4) m.fs.unit.control_volume.properties_out[0.0].pressure_sat_comp.setub(5e6) + m.fs.properties.set_default_scaling("flow_mol", 1) + m.fs.properties.set_default_scaling("flow_mol_phase", 1) + m.fs.properties.set_default_scaling("flow_mol_phase_comp", 1) + + return m @pytest.mark.build @@ -291,7 +297,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -496,7 +507,12 @@ def test_conservation(self, iapws): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory('core.scale_model').create_using( + iapws, + rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_heat_exchanger.py b/idaes/models/unit_models/tests/test_heat_exchanger.py index 85f1797fb2..0b8887d37a 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger.py @@ -299,6 +299,9 @@ def basic_model(cb=delta_temperature_lmtd_callback): m.fs.unit.area.fix(1000) m.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(m.fs.unit.area, 1e-3) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + assert degrees_of_freedom(m) == 0 m.fs.unit.initialize() return m @@ -328,6 +331,9 @@ def basic_model2(cb=delta_temperature_lmtd_callback): m.fs.unit.area.fix(100) m.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(m.fs.unit.area, 1e-2) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + assert degrees_of_freedom(m) == 0 m.fs.unit.initialize() return m @@ -357,6 +363,10 @@ def basic_model3(cb=delta_temperature_lmtd_callback): m.fs.unit.area.fix(1000) m.fs.unit.overall_heat_transfer_coefficient.fix(100) m.fs.unit.crossflow_factor.fix(1.0) + + iscale.set_scaling_factor(m.fs.unit.area, 1e-3) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + assert degrees_of_freedom(m) == 0 m.fs.unit.initialize() return m @@ -503,6 +513,9 @@ def btx(self): m.fs.unit.area.fix(1) m.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(m.fs.unit.area, 1) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + # Bound temperature differences to avoid division by zero m.fs.unit.delta_temperature_in[0.0].setlb(40) m.fs.unit.delta_temperature_out[0.0].setlb(0.1) @@ -714,7 +727,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -750,6 +768,9 @@ def btx(self): m.fs.unit.area.fix(1) m.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(m.fs.unit.area, 1) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + # Bound temperature differences to avoid division by zero m.fs.unit.delta_temperature_in[0.0].setlb(40) m.fs.unit.delta_temperature_out[0.0].setlb(0.1) @@ -956,7 +977,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -988,6 +1014,9 @@ def iapws(self): m.fs.unit.area.fix(1000) m.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(m.fs.unit.area, 1e-3) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + return m @pytest.fixture(scope="class") @@ -1015,6 +1044,9 @@ def iapws_underwood(self): m.fs.unit.area.fix(1000) m.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(m.fs.unit.area, 1e-3) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + return m @pytest.mark.build @@ -1075,6 +1107,9 @@ def test_dof_alt_name1(self, iapws): iapws.fs.unit.area.fix(1000) iapws.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(iapws.fs.unit.area, 1e-3) + iscale.set_scaling_factor(iapws.fs.unit.overall_heat_transfer_coefficient, 1e-2) + assert degrees_of_freedom(iapws) == 0 @pytest.mark.ui @@ -1246,10 +1281,10 @@ def test_conservation(self, iapws): @pytest.mark.component def test_numerical_issues(self, iapws): iscale.calculate_scaling_factors(iapws) - m_scaled = TransformationFactory("core.scale_model").create_using( + iapws_scaled = TransformationFactory("core.scale_model").create_using( iapws, rename=False ) - dt = DiagnosticsToolbox(m_scaled) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() @@ -1290,6 +1325,9 @@ def sapon(self): m.fs.unit.overall_heat_transfer_coefficient.fix(100) m.fs.unit.crossflow_factor.fix(0.6) + iscale.set_scaling_factor(m.fs.unit.area, 1e-3) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + return m @pytest.mark.build @@ -1535,7 +1573,7 @@ def test_numerical_issues(self, sapon): # Heat transfer constraint has a a residual of ~1e-3 # Model could be better scaled dt = DiagnosticsToolbox(sapon, constraint_residual_tolerance=1e-2) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) # ----------------------------------------------------------------------------- @@ -1688,6 +1726,9 @@ def btx(self): m.fs.unit.cold_side.scaling_factor_pressure = 1 + iscale.set_scaling_factor(m.fs.unit.area, 1) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + return m @pytest.mark.build @@ -1898,7 +1939,13 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) + dt.assert_no_numerical_warnings() @pytest.mark.component @@ -2055,6 +2102,9 @@ def model(self): m.fs.unit.area.fix(1) m.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(m.fs.unit.area, 1) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + m.fs.unit.cold_side.scaling_factor_pressure = 1 return m @@ -2139,6 +2189,9 @@ def model(self): m.fs.unit.area.fix(1000) m.fs.unit.overall_heat_transfer_coefficient.fix(100) + iscale.set_scaling_factor(m.fs.unit.area, 1e-3) + iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) + return m @pytest.mark.component diff --git a/idaes/models/unit_models/tests/test_heat_exchanger_1D.py b/idaes/models/unit_models/tests/test_heat_exchanger_1D.py index 2724dc6fc3..b6e4076295 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger_1D.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger_1D.py @@ -732,7 +732,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -964,7 +969,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -1562,7 +1572,11 @@ def test_conservation(self, iapws): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() @@ -1779,7 +1793,11 @@ def test_conservation(self, iapws): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() @@ -2038,7 +2056,7 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) # # ----------------------------------------------------------------------------- @@ -2296,7 +2314,7 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) # # ----------------------------------------------------------------------------- @@ -2640,7 +2658,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.integration def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @pytest.mark.component From 330953fcc6016a76882c150c86fd880064131a63 Mon Sep 17 00:00:00 2001 From: Doug A Date: Mon, 17 Jun 2024 14:01:06 -0400 Subject: [PATCH 15/28] finish cleanup --- idaes/core/base/control_volume1d.py | 26 ++++++++++++++ .../activity_coeff_prop_pack.py | 36 +++++++++++++------ idaes/models/unit_models/separator.py | 28 +++------------ .../tests/test_equilibrium_reactor.py | 2 +- idaes/models/unit_models/tests/test_heater.py | 26 +++++++++++--- idaes/models/unit_models/tests/test_hx_ntu.py | 9 ++++- idaes/models/unit_models/tests/test_mixer.py | 11 ++++-- idaes/models/unit_models/tests/test_pfr.py | 2 +- .../tests/test_pressure_changer.py | 16 +++++++-- .../unit_models/tests/test_separator.py | 25 ++++++++++--- .../tests/test_shell_and_tube_1D.py | 35 +++++++++++++++--- .../column_models/tests/test_reboiler.py | 2 +- 12 files changed, 161 insertions(+), 57 deletions(-) diff --git a/idaes/core/base/control_volume1d.py b/idaes/core/base/control_volume1d.py index 733d8943dd..045415b71e 100644 --- a/idaes/core/base/control_volume1d.py +++ b/idaes/core/base/control_volume1d.py @@ -2554,3 +2554,29 @@ 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) diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index 28a6eea414..bf73c3be5e 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -1562,30 +1562,30 @@ def get_enthalpy_flow_terms(self, p): def rule_enthalpy_flow_terms(b, p): if p == "Vap": - if self.params.config.state_vars == "FTPz": + if b.params.config.state_vars == "FTPz": return ( - self.flow_mol_phase["Vap"] * self.enth_mol_phase["Vap"] + b.flow_mol_phase["Vap"] * b.enth_mol_phase["Vap"] ) else: return ( sum( - self.flow_mol_phase_comp["Vap", j] - for j in self.params.component_list + b.flow_mol_phase_comp["Vap", j] + for j in b.params.component_list ) - * self.enth_mol_phase["Vap"] + * b.enth_mol_phase["Vap"] ) elif p == "Liq": - if self.params.config.state_vars == "FTPz": + if b.params.config.state_vars == "FTPz": return ( - self.flow_mol_phase["Liq"] * self.enth_mol_phase["Liq"] + b.flow_mol_phase["Liq"] * b.enth_mol_phase["Liq"] ) else: return ( sum( - self.flow_mol_phase_comp["Liq", j] - for j in self.params.component_list + b.flow_mol_phase_comp["Liq", j] + for j in b.params.component_list ) - * self.enth_mol_phase["Liq"] + * b.enth_mol_phase["Liq"] ) self.enthalpy_flow_terms = Expression( @@ -2238,3 +2238,19 @@ def calculate_scaling_factors(self): self.eq_P_vap[p], default=1, warning=True ) iscale.constraint_scaling_transform(c, sf, overwrite=False) + + if self.is_property_constructed("enthalpy_flow_terms"): + for p, expr in self.enthalpy_flow_terms.items(): + sf_enth_mol = iscale.get_scaling_factor(self.enth_mol_phase[p], default=1) + if self.params.config.state_vars == "FTPz": + sf_flow = iscale.get_scaling_factor(self.flow_mol_phase[p], default=1) + else: + inv_sf_flow = 0 + for j in self.params.component_list: + inv_sf_flow += 1/iscale.get_scaling_factor( + self.flow_mol_phase_comp[p,j], + default=1 + ) + sf_flow = 1/inv_sf_flow + + iscale.set_scaling_factor(expr, sf_enth_mol*sf_flow, overwrite=False) diff --git a/idaes/models/unit_models/separator.py b/idaes/models/unit_models/separator.py index 168f30aded..9bf1c503a0 100644 --- a/idaes/models/unit_models/separator.py +++ b/idaes/models/unit_models/separator.py @@ -1829,33 +1829,21 @@ def calculate_scaling_factors(self): s = _s if _s < s else s iscale.constraint_scaling_transform(c, s, overwrite=False) - if self.config.energy_split_basis == EnergySplittingType.none: - pass - elif self.config.energy_split_basis == EnergySplittingType.equal_temperature: + if hasattr(self, "temperature_equality_eqn"): for (t, i), c in self.temperature_equality_eqn.items(): s = iscale.get_scaling_factor( mixed_state[t].temperature, default=1, warning=True ) iscale.constraint_scaling_transform(c, s, overwrite=False) - elif self.config.energy_split_basis == EnergySplittingType.equal_molar_enthalpy: - for (t, i), c in self.temperature_equality_eqn.items(): + if hasattr(self, "molar_enthalpy_equality_eqn"): + for (t, i), c in self.molar_enthalpy_equality_eqn.items(): s_enth_mol = iscale.get_scaling_factor( mixed_state[t].enth_mol, default=1, warning=True ) iscale.constraint_scaling_transform(c, s_enth_mol, overwrite=False) - elif self.config.energy_split_basis == EnergySplittingType.enthalpy_split: - # Validate split fraction type - if ( - self.config.split_basis == SplittingType.phaseComponentFlow - or self.config.split_basis == SplittingType.componentFlow - ): - raise ConfigurationError( - "{} Cannot use energy_split_basis == enthalpy_split " - "with split_basis == component or phaseComponent.".format(self.name) - ) - + if hasattr(self, "molar_enthalpy_splitting_eqn"): for (t, i), c in self.molar_enthalpy_splitting_eqn.items(): sf_enth = float('inf') for p in mixed_state[t].phase_list: @@ -1869,14 +1857,6 @@ def calculate_scaling_factors(self): ) iscale.constraint_scaling_transform(c, sf_enth, overwrite=False) - else: - raise BurntToast( - "{} received unrecognised value for " - "energy_split_basis. This should never happen, so" - " please contact the IDAES developers with this " - "bug.".format(self.name) - ) - def _get_performance_contents(self, time_point=0): if hasattr(self, "split_fraction"): var_dict = {} diff --git a/idaes/models/unit_models/tests/test_equilibrium_reactor.py b/idaes/models/unit_models/tests/test_equilibrium_reactor.py index 8dc8c9c199..8a90cafe00 100644 --- a/idaes/models/unit_models/tests/test_equilibrium_reactor.py +++ b/idaes/models/unit_models/tests/test_equilibrium_reactor.py @@ -238,7 +238,7 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.ui @pytest.mark.unit diff --git a/idaes/models/unit_models/tests/test_heater.py b/idaes/models/unit_models/tests/test_heater.py index dbee5e970b..20ea6cca1c 100644 --- a/idaes/models/unit_models/tests/test_heater.py +++ b/idaes/models/unit_models/tests/test_heater.py @@ -20,6 +20,7 @@ from pyomo.environ import ( check_optimal_termination, ConcreteModel, + TransformationFactory, value, units as pyunits, ) @@ -58,6 +59,7 @@ InitializationStatus, ) from idaes.core.util import DiagnosticsToolbox +import idaes.core.util.scaling as iscale # ----------------------------------------------------------------------------- # Get default solver for testing @@ -188,7 +190,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @pytest.mark.ui @@ -218,6 +225,8 @@ def iapws(self): m.fs.unit.heat_duty.fix(10000) + m.fs.properties.set_default_scaling("flow_mol", 1) + return m @pytest.mark.build @@ -299,7 +308,11 @@ def test_conservation(self, iapws): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() @pytest.mark.ui @@ -480,7 +493,7 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.ui @pytest.mark.unit @@ -594,7 +607,12 @@ def test_conservation(self, btg): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btg): - dt = DiagnosticsToolbox(btg) + iscale.calculate_scaling_factors(btg) + btg_scaled = TransformationFactory('core.scale_model').create_using( + btg, + rename=False + ) + dt = DiagnosticsToolbox(btg_scaled) dt.assert_no_numerical_warnings() @pytest.mark.ui diff --git a/idaes/models/unit_models/tests/test_hx_ntu.py b/idaes/models/unit_models/tests/test_hx_ntu.py index 86af622fa1..14fc0e98dc 100644 --- a/idaes/models/unit_models/tests/test_hx_ntu.py +++ b/idaes/models/unit_models/tests/test_hx_ntu.py @@ -23,6 +23,7 @@ Constraint, Expression, Param, + TransformationFactory, units as pyunits, value, Var, @@ -50,6 +51,7 @@ InitializationStatus, ) from idaes.core.util import DiagnosticsToolbox +import idaes.core.util.scaling as iscale # ----------------------------------------------------------------------------- @@ -501,7 +503,12 @@ def test_conservation(self, model): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, model): - dt = DiagnosticsToolbox(model) + iscale.calculate_scaling_factors(model) + model_scaled = TransformationFactory('core.scale_model').create_using( + model, + rename=False + ) + dt = DiagnosticsToolbox(model_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_mixer.py b/idaes/models/unit_models/tests/test_mixer.py index 6046bf23e8..f5c0e0c100 100644 --- a/idaes/models/unit_models/tests/test_mixer.py +++ b/idaes/models/unit_models/tests/test_mixer.py @@ -25,6 +25,7 @@ Param, RangeSet, Set, + TransformationFactory, Var, value, units as pyunits, @@ -90,6 +91,7 @@ fix_state_vars, ) from idaes.core.util import DiagnosticsToolbox +import idaes.core.util.scaling as iscale # TODO: Should have a test for this that does not require models_extra @@ -936,7 +938,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -1462,7 +1469,7 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.skipif(not cubic_roots_available(), reason="Cubic functions not available") diff --git a/idaes/models/unit_models/tests/test_pfr.py b/idaes/models/unit_models/tests/test_pfr.py index b639ba0650..faf4fbf60c 100644 --- a/idaes/models/unit_models/tests/test_pfr.py +++ b/idaes/models/unit_models/tests/test_pfr.py @@ -263,7 +263,7 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.ui @pytest.mark.unit diff --git a/idaes/models/unit_models/tests/test_pressure_changer.py b/idaes/models/unit_models/tests/test_pressure_changer.py index 0d2ab161d6..4bf73eb1a5 100644 --- a/idaes/models/unit_models/tests/test_pressure_changer.py +++ b/idaes/models/unit_models/tests/test_pressure_changer.py @@ -21,6 +21,7 @@ check_optimal_termination, ConcreteModel, Constraint, + TransformationFactory, units, value, Var, @@ -337,7 +338,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @pytest.mark.ui @@ -571,7 +577,11 @@ def test_verify(self, iapws_turb): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() @pytest.mark.ui @@ -734,7 +744,7 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.ui @pytest.mark.unit diff --git a/idaes/models/unit_models/tests/test_separator.py b/idaes/models/unit_models/tests/test_separator.py index 13f9458448..adb24c57b8 100644 --- a/idaes/models/unit_models/tests/test_separator.py +++ b/idaes/models/unit_models/tests/test_separator.py @@ -24,6 +24,7 @@ Constraint, Param, Set, + TransformationFactory, value, Var, units as pyunits, @@ -1114,7 +1115,7 @@ def test_conservation(self, sapon): @pytest.mark.component def test_numerical_issues(self, sapon): dt = DiagnosticsToolbox(sapon) - dt.assert_no_numerical_warnings() + dt.assert_no_numerical_warnings(ignore_parallel_components=True) # ----------------------------------------------------------------------------- @@ -1419,7 +1420,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -1632,7 +1638,11 @@ def test_conservation(self, iapws): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() @@ -3252,8 +3262,13 @@ def test_conservation(self, btx): @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component - def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + def test_numerical_issues(self, btx): + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_shell_and_tube_1D.py b/idaes/models/unit_models/tests/test_shell_and_tube_1D.py index 362e267313..072936806d 100644 --- a/idaes/models/unit_models/tests/test_shell_and_tube_1D.py +++ b/idaes/models/unit_models/tests/test_shell_and_tube_1D.py @@ -20,6 +20,7 @@ from pyomo.environ import ( check_optimal_termination, ConcreteModel, + TransformationFactory, value, units as pyunits, ) @@ -548,10 +549,16 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() + # ----------------------------------------------------------------------------- class TestBTX_countercurrent(object): @pytest.fixture(scope="class") @@ -793,7 +800,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -1028,7 +1040,11 @@ def test_conservation(self, iapws): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() @@ -1263,7 +1279,11 @@ def test_conservation(self, iapws): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component def test_numerical_issues(self, iapws): - dt = DiagnosticsToolbox(iapws) + iscale.calculate_scaling_factors(iapws) + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False + ) + dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() @@ -1624,7 +1644,12 @@ def test_conservation(self, btx): @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.integration def test_numerical_issues(self, btx): - dt = DiagnosticsToolbox(btx) + iscale.calculate_scaling_factors(btx) + btx_scaled = TransformationFactory('core.scale_model').create_using( + btx, + rename=False + ) + dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @pytest.mark.component diff --git a/idaes/models_extra/column_models/tests/test_reboiler.py b/idaes/models_extra/column_models/tests/test_reboiler.py index fac8567fcb..62c44c41da 100644 --- a/idaes/models_extra/column_models/tests/test_reboiler.py +++ b/idaes/models_extra/column_models/tests/test_reboiler.py @@ -281,7 +281,7 @@ def test_solution(self, btx_ftpz, btx_fctp): ) # Unit level - assert pytest.approx(16926.5, rel=1e-4) == value(btx_fctp.fs.unit.heat_duty[0]) + assert pytest.approx(16916.9, rel=1e-4) == value(btx_fctp.fs.unit.heat_duty[0]) @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component From e30c6d5f884c2f628618dff3e247d144ca051aef Mon Sep 17 00:00:00 2001 From: Doug A Date: Mon, 17 Jun 2024 14:41:12 -0400 Subject: [PATCH 16/28] commit black --- idaes/core/base/control_volume1d.py | 6 ++-- .../activity_coeff_prop_pack.py | 35 +++++++++---------- idaes/models/unit_models/feed_flash.py | 19 +++++----- idaes/models/unit_models/separator.py | 6 ++-- .../unit_models/tests/test_feed_flash.py | 12 +++---- idaes/models/unit_models/tests/test_flash.py | 13 +++---- .../unit_models/tests/test_heat_exchanger.py | 15 ++++---- .../tests/test_heat_exchanger_1D.py | 15 ++++---- idaes/models/unit_models/tests/test_heater.py | 10 +++--- idaes/models/unit_models/tests/test_hx_ntu.py | 5 ++- idaes/models/unit_models/tests/test_mixer.py | 5 ++- .../tests/test_pressure_changer.py | 5 ++- .../unit_models/tests/test_separator.py | 12 +++---- .../tests/test_shell_and_tube_1D.py | 16 ++++----- 14 files changed, 75 insertions(+), 99 deletions(-) diff --git a/idaes/core/base/control_volume1d.py b/idaes/core/base/control_volume1d.py index 045415b71e..265c36907c 100644 --- a/idaes/core/base/control_volume1d.py +++ b/idaes/core/base/control_volume1d.py @@ -2569,12 +2569,10 @@ def calculate_scaling_factors(self): 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 + 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 diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index bf73c3be5e..06e096da9a 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -222,7 +222,7 @@ def build(self): self.set_default_scaling("temperature_dew", 1e-1) self.set_default_scaling("temperature_bubble", 1e-1) self.set_default_scaling("flow_mol_phase", 1e-2) - self.set_default_scaling("density_mol", 1/11.1e3, index="Liq") + self.set_default_scaling("density_mol", 1 / 11.1e3, index="Liq") self.set_default_scaling("density_mol", 0.31, index="Vap") self.set_default_scaling("pressure", 1e-6) @@ -236,8 +236,8 @@ def build(self): self.set_default_scaling("enth_mol_phase", 1e-4, index="Vap") self.set_default_scaling("entr_mol_phase_comp", 1e-2) self.set_default_scaling("entr_mol_phase", 1e-2) - self.set_default_scaling("mw", 1/100) - self.set_default_scaling("mw_phase", 1/100) + self.set_default_scaling("mw", 1 / 100) + self.set_default_scaling("mw_phase", 1 / 100) self.set_default_scaling("gibbs_mol_phase_comp", 1e-3) self.set_default_scaling("fug_vap", 1e-6) self.set_default_scaling("fug_liq", 1e-6) @@ -1563,9 +1563,7 @@ def get_enthalpy_flow_terms(self, p): def rule_enthalpy_flow_terms(b, p): if p == "Vap": if b.params.config.state_vars == "FTPz": - return ( - b.flow_mol_phase["Vap"] * b.enth_mol_phase["Vap"] - ) + return b.flow_mol_phase["Vap"] * b.enth_mol_phase["Vap"] else: return ( sum( @@ -1576,9 +1574,7 @@ def rule_enthalpy_flow_terms(b, p): ) elif p == "Liq": if b.params.config.state_vars == "FTPz": - return ( - b.flow_mol_phase["Liq"] * b.enth_mol_phase["Liq"] - ) + return b.flow_mol_phase["Liq"] * b.enth_mol_phase["Liq"] else: return ( sum( @@ -2230,7 +2226,7 @@ def calculate_scaling_factors(self): for j in self.mole_frac_comp } for j, c in self.eq_phase_equilibrium.items(): - iscale.constraint_scaling_transform(c, sf_p*sf_x[j], overwrite=False) + iscale.constraint_scaling_transform(c, sf_p * sf_x[j], overwrite=False) if self.is_property_constructed("eq_P_vap"): for p, c in self.eq_P_vap.items(): @@ -2240,17 +2236,20 @@ def calculate_scaling_factors(self): iscale.constraint_scaling_transform(c, sf, overwrite=False) if self.is_property_constructed("enthalpy_flow_terms"): - for p, expr in self.enthalpy_flow_terms.items(): - sf_enth_mol = iscale.get_scaling_factor(self.enth_mol_phase[p], default=1) + for p, expr in self.enthalpy_flow_terms.items(): + sf_enth_mol = iscale.get_scaling_factor( + self.enth_mol_phase[p], default=1 + ) if self.params.config.state_vars == "FTPz": - sf_flow = iscale.get_scaling_factor(self.flow_mol_phase[p], default=1) + sf_flow = iscale.get_scaling_factor( + self.flow_mol_phase[p], default=1 + ) else: inv_sf_flow = 0 for j in self.params.component_list: - inv_sf_flow += 1/iscale.get_scaling_factor( - self.flow_mol_phase_comp[p,j], - default=1 + inv_sf_flow += 1 / iscale.get_scaling_factor( + self.flow_mol_phase_comp[p, j], default=1 ) - sf_flow = 1/inv_sf_flow + sf_flow = 1 / inv_sf_flow - iscale.set_scaling_factor(expr, sf_enth_mol*sf_flow, overwrite=False) + iscale.set_scaling_factor(expr, sf_enth_mol * sf_flow, overwrite=False) diff --git a/idaes/models/unit_models/feed_flash.py b/idaes/models/unit_models/feed_flash.py index 1cec7fefc7..7bde834f7c 100644 --- a/idaes/models/unit_models/feed_flash.py +++ b/idaes/models/unit_models/feed_flash.py @@ -223,9 +223,11 @@ def calculate_scaling_factors(self): sT = iscale.get_scaling_factor( self.control_volume.properties_in[t].temperature, default=1, - warning=True + warning=True, + ) + iscale.constraint_scaling_transform( + self.isothermal[t], sT, overwrite=False ) - iscale.constraint_scaling_transform(self.isothermal[t], sT, overwrite=False) elif self.config.flash_type == FlashType.isenthalpic: cv = self.control_volume for t in self.flowsheet().time: @@ -235,11 +237,10 @@ def calculate_scaling_factors(self): s_enth, iscale.get_scaling_factor( cv.properties_in[t].get_enthalpy_flow_terms(p), - default = 1, - warning=True - ) + default=1, + warning=True, + ), ) - iscale.constraint_scaling_transform(self.isenthalpic[t], s_enth, overwrite=False) - - - + iscale.constraint_scaling_transform( + self.isenthalpic[t], s_enth, overwrite=False + ) diff --git a/idaes/models/unit_models/separator.py b/idaes/models/unit_models/separator.py index 9bf1c503a0..7e97d7fce0 100644 --- a/idaes/models/unit_models/separator.py +++ b/idaes/models/unit_models/separator.py @@ -1845,15 +1845,15 @@ def calculate_scaling_factors(self): if hasattr(self, "molar_enthalpy_splitting_eqn"): for (t, i), c in self.molar_enthalpy_splitting_eqn.items(): - sf_enth = float('inf') + sf_enth = float("inf") for p in mixed_state[t].phase_list: sf_enth = min( sf_enth, iscale.get_scaling_factor( mixed_state[t].get_enthalpy_flow_terms(p), default=1, - warning=True - ) + warning=True, + ), ) iscale.constraint_scaling_transform(c, sf_enth, overwrite=False) diff --git a/idaes/models/unit_models/tests/test_feed_flash.py b/idaes/models/unit_models/tests/test_feed_flash.py index 23aa2b4fd1..302cea5280 100644 --- a/idaes/models/unit_models/tests/test_feed_flash.py +++ b/idaes/models/unit_models/tests/test_feed_flash.py @@ -216,9 +216,8 @@ def test_solution(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -340,11 +339,10 @@ def test_solution(self, iapws): @pytest.mark.component def test_numerical_issues(self, iapws): dt = DiagnosticsToolbox(iapws) - + iscale.calculate_scaling_factors(iapws) - iapws_scaled = TransformationFactory('core.scale_model').create_using( - iapws, - rename=False + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False ) dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_flash.py b/idaes/models/unit_models/tests/test_flash.py index 02e621635f..490e5406f6 100644 --- a/idaes/models/unit_models/tests/test_flash.py +++ b/idaes/models/unit_models/tests/test_flash.py @@ -102,7 +102,7 @@ def btx(self): m.fs.properties = BTXParameterBlock( valid_phase=("Liq", "Vap"), activity_coeff_model="Ideal" ) - + m.fs.unit = Flash(property_package=m.fs.properties) m.fs.unit.inlet.flow_mol.fix(1) @@ -139,7 +139,6 @@ def btx(self): m.fs.properties.set_default_scaling("flow_mol_phase", 1) m.fs.properties.set_default_scaling("flow_mol_phase_comp", 1) - return m @pytest.mark.build @@ -298,9 +297,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -508,9 +506,8 @@ def test_conservation(self, iapws): @pytest.mark.component def test_numerical_issues(self, iapws): iscale.calculate_scaling_factors(iapws) - iapws_scaled = TransformationFactory('core.scale_model').create_using( - iapws, - rename=False + iapws_scaled = TransformationFactory("core.scale_model").create_using( + iapws, rename=False ) dt = DiagnosticsToolbox(iapws_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_heat_exchanger.py b/idaes/models/unit_models/tests/test_heat_exchanger.py index 0b8887d37a..a41440a42f 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger.py @@ -728,9 +728,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -978,9 +977,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -1940,9 +1938,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) diff --git a/idaes/models/unit_models/tests/test_heat_exchanger_1D.py b/idaes/models/unit_models/tests/test_heat_exchanger_1D.py index b6e4076295..a6157de853 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger_1D.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger_1D.py @@ -733,9 +733,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -970,9 +969,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -2659,9 +2657,8 @@ def test_conservation(self, btx): @pytest.mark.integration def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_heater.py b/idaes/models/unit_models/tests/test_heater.py index 20ea6cca1c..aaf9e4881a 100644 --- a/idaes/models/unit_models/tests/test_heater.py +++ b/idaes/models/unit_models/tests/test_heater.py @@ -191,9 +191,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -608,9 +607,8 @@ def test_conservation(self, btg): @pytest.mark.component def test_numerical_issues(self, btg): iscale.calculate_scaling_factors(btg) - btg_scaled = TransformationFactory('core.scale_model').create_using( - btg, - rename=False + btg_scaled = TransformationFactory("core.scale_model").create_using( + btg, rename=False ) dt = DiagnosticsToolbox(btg_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_hx_ntu.py b/idaes/models/unit_models/tests/test_hx_ntu.py index 14fc0e98dc..ecf06551f0 100644 --- a/idaes/models/unit_models/tests/test_hx_ntu.py +++ b/idaes/models/unit_models/tests/test_hx_ntu.py @@ -504,9 +504,8 @@ def test_conservation(self, model): @pytest.mark.component def test_numerical_issues(self, model): iscale.calculate_scaling_factors(model) - model_scaled = TransformationFactory('core.scale_model').create_using( - model, - rename=False + model_scaled = TransformationFactory("core.scale_model").create_using( + model, rename=False ) dt = DiagnosticsToolbox(model_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_mixer.py b/idaes/models/unit_models/tests/test_mixer.py index f5c0e0c100..3132758b26 100644 --- a/idaes/models/unit_models/tests/test_mixer.py +++ b/idaes/models/unit_models/tests/test_mixer.py @@ -939,9 +939,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_pressure_changer.py b/idaes/models/unit_models/tests/test_pressure_changer.py index 4bf73eb1a5..42c13acffa 100644 --- a/idaes/models/unit_models/tests/test_pressure_changer.py +++ b/idaes/models/unit_models/tests/test_pressure_changer.py @@ -339,9 +339,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_separator.py b/idaes/models/unit_models/tests/test_separator.py index adb24c57b8..2d21f98953 100644 --- a/idaes/models/unit_models/tests/test_separator.py +++ b/idaes/models/unit_models/tests/test_separator.py @@ -1421,9 +1421,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -3262,11 +3261,10 @@ def test_conservation(self, btx): @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component - def test_numerical_issues(self, btx): + def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() diff --git a/idaes/models/unit_models/tests/test_shell_and_tube_1D.py b/idaes/models/unit_models/tests/test_shell_and_tube_1D.py index 072936806d..4012191b84 100644 --- a/idaes/models/unit_models/tests/test_shell_and_tube_1D.py +++ b/idaes/models/unit_models/tests/test_shell_and_tube_1D.py @@ -550,15 +550,13 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() - # ----------------------------------------------------------------------------- class TestBTX_countercurrent(object): @pytest.fixture(scope="class") @@ -801,9 +799,8 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() @@ -1645,9 +1642,8 @@ def test_conservation(self, btx): @pytest.mark.integration def test_numerical_issues(self, btx): iscale.calculate_scaling_factors(btx) - btx_scaled = TransformationFactory('core.scale_model').create_using( - btx, - rename=False + btx_scaled = TransformationFactory("core.scale_model").create_using( + btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) dt.assert_no_numerical_warnings() From 50903cc7164e2688d6e7be87a0626ee633e24bea Mon Sep 17 00:00:00 2001 From: Doug A Date: Mon, 17 Jun 2024 15:01:22 -0400 Subject: [PATCH 17/28] pylint --- idaes/models/unit_models/feed_flash.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/idaes/models/unit_models/feed_flash.py b/idaes/models/unit_models/feed_flash.py index 7bde834f7c..c8b84095af 100644 --- a/idaes/models/unit_models/feed_flash.py +++ b/idaes/models/unit_models/feed_flash.py @@ -16,7 +16,7 @@ from enum import Enum # Import Pyomo libraries -from pyomo.environ import exp, log, Reference +from pyomo.environ import Reference from pyomo.common.config import ConfigBlock, ConfigValue, In # Import IDAES cores From a54f189048af66d59336dbfa96c2400ea7e3f995 Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 18 Jun 2024 09:56:45 -0400 Subject: [PATCH 18/28] Begun dealing with negative phase flows --- idaes/models/flowsheets/demo_flowsheet.py | 391 +----------------- .../flowsheets/tests/test_demo_flowsheet.py | 2 + .../activity_coeff_prop_pack.py | 21 +- 3 files changed, 36 insertions(+), 378 deletions(-) diff --git a/idaes/models/flowsheets/demo_flowsheet.py b/idaes/models/flowsheets/demo_flowsheet.py index 7708ae442b..6cc73f9e13 100644 --- a/idaes/models/flowsheets/demo_flowsheet.py +++ b/idaes/models/flowsheets/demo_flowsheet.py @@ -56,15 +56,22 @@ def build_flowsheet(): def set_scaling(m): """Set scaling for demo flowsheet""" - - iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].flow_mol_phase, 1e6) + params = m.fs.BT_props + params.set_default_scaling("flow_mol", 1) + params.set_default_scaling("flow_mol_phase", 1) + params.set_default_scaling("flow_mol_phase_comp", 1) + params.set_default_scaling("mole_frac_comp", 1) + params.set_default_scaling("mole_frac_phase_comp", 1) + params.set_default_scaling("temperature", 1e-1) + params.set_default_scaling("pressure", 1e-5) + params.set_default_scaling("pressure_sat_comp", 1e-5, index="benzene") + params.set_default_scaling("pressure_sat_comp", 1e-4, index="toluene") + params.set_default_scaling("enth_mol_phase_comp", 1e-4) + + # Mixer M01 iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].flow_mol_phase["Liq"], 1e6) - iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].flow_mol_phase["Vap"], 1) - iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].mole_frac_phase_comp["Liq", "benzene"], 1 - ) iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].mole_frac_phase_comp["Vap", "benzene"], 1 + m.fs.M01.inlet_1_state[0].mole_frac_comp["toluene"], 1e5 ) iscale.set_scaling_factor( m.fs.M01.inlet_1_state[0].mole_frac_phase_comp["Liq", "toluene"], 1e5 @@ -73,380 +80,20 @@ def set_scaling(m): m.fs.M01.inlet_1_state[0].mole_frac_phase_comp["Vap", "toluene"], 1e6 ) iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].pressure_sat_comp["benzene"], 1e-5 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].pressure_sat_comp["toluene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].enth_mol_phase_comp["Liq", "benzene"], 1e-4 + m.fs.M01.inlet_2_state[0].mole_frac_comp["benzene"], 1e5 ) - iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].enth_mol_phase_comp["Vap", "benzene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].enth_mol_phase_comp["Liq", "toluene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].enth_mol_phase_comp["Vap", "toluene"], 1e-4 - ) - - iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].flow_mol_phase, 1e1) - iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].flow_mol_phase["Vap"], 1e1) - iscale.set_scaling_factor( m.fs.M01.inlet_2_state[0].mole_frac_phase_comp["Liq", "benzene"], 1e5 ) iscale.set_scaling_factor( m.fs.M01.inlet_2_state[0].mole_frac_phase_comp["Vap", "benzene"], 1e5 ) - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].mole_frac_phase_comp["Liq", "toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].mole_frac_phase_comp["Vap", "toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].pressure_sat_comp["benzene"], 1e-6 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].pressure_sat_comp["toluene"], 1e-6 - ) - - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].enth_mol_phase_comp["Liq", "benzene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].enth_mol_phase_comp["Vap", "benzene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].enth_mol_phase_comp["Liq", "toluene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].enth_mol_phase_comp["Vap", "toluene"], 1e-4 - ) + # Heater H02 + iscale.set_scaling_factor(m.fs.H02.control_volume.heat[0], 1e-5) - iscale.set_scaling_factor( - m.fs.M01.mixed_state[0].enth_mol_phase_comp["Liq", "benzene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.mixed_state[0].enth_mol_phase_comp["Vap", "benzene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.mixed_state[0].enth_mol_phase_comp["Liq", "toluene"], 1e-4 - ) - iscale.set_scaling_factor( - m.fs.M01.mixed_state[0].enth_mol_phase_comp["Vap", "toluene"], 1e-4 - ) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].enthalpy_flow_terms["Liq"], 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].enthalpy_flow_terms["Vap"], 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_total, 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_sum_mol_frac, 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_mol_frac_out, 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_comp["benzene"], 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_comp["toluene"], 1) - iscale.set_scaling_factor( - m.fs.M01.mixed_state[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.M01.mixed_state[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_P_vap["benzene"], 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].eq_P_vap["toluene"], 1) - - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].pressure, 1e-5) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].pressure_sat_comp, 1e-5) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].temperature, 1e-2) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].flow_mol, 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].flow_mol_phase, 1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].mole_frac_comp["benzene"], 1e1) - iscale.set_scaling_factor(m.fs.M01.mixed_state[0].mole_frac_comp["toluene"], 1e1) - - iscale.set_scaling_factor(m.fs.H02.control_volume.heat, 1e-5) - iscale.set_scaling_factor(m.fs.H02.control_volume.properties_in[0].pressure, 1e-5) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].pressure_sat_comp, 1e-5 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].temperature, 1e-2 - ) - iscale.set_scaling_factor(m.fs.H02.control_volume.properties_in[0].flow_mol, 1) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].flow_mol_phase, 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].mole_frac_comp["benzene"], 1e1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].mole_frac_comp["toluene"], 1e1 - ) - iscale.set_scaling_factor(m.fs.H02.control_volume.properties_out[0].flow_mol, 1) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].flow_mol_phase, 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].mole_frac_comp["benzene"], 1e1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].mole_frac_comp["toluene"], 1e1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].temperature, 1e-2 - ) - iscale.set_scaling_factor(m.fs.H02.control_volume.properties_out[0].pressure, 1e-5) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].pressure_sat_comp, 1e-5 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].enth_mol_phase_comp["Liq", "benzene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].enth_mol_phase_comp["Vap", "benzene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].enth_mol_phase_comp["Liq", "toluene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].enth_mol_phase_comp["Vap", "toluene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].eq_comp["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].eq_comp["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].eq_comp["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].eq_comp["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].eq_P_vap["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].eq_P_vap["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0.0].eq_P_vap["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0.0].eq_P_vap["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].enth_mol_phase_comp["Liq", "benzene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].enth_mol_phase_comp["Vap", "benzene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].enth_mol_phase_comp["Liq", "toluene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0].enth_mol_phase_comp["Vap", "toluene"], - 1e-4, - ) - iscale.set_scaling_factor(m.fs.H02.control_volume.properties_in[0.0].eq_total, 1) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].eq_sum_mol_frac, 1 - ) - iscale.set_scaling_factor(m.fs.H02.control_volume.properties_out[0.0].eq_total, 1) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0.0].eq_sum_mol_frac, 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_out[0.0].eq_mol_frac_out, 1 - ) - - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].material_flow_terms[ - "Liq", "benzene" - ], - 1, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].material_flow_terms[ - "Vap", "benzene" - ], - 1, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].material_flow_terms[ - "Liq", "toluene" - ], - 1, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].material_flow_terms[ - "Vap", "toluene" - ], - 1, - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].enthalpy_flow_terms["Liq"], 1 - ) - iscale.set_scaling_factor( - m.fs.H02.control_volume.properties_in[0.0].enthalpy_flow_terms["Vap"], 1 - ) - - iscale.set_scaling_factor(m.fs.F03.control_volume.heat, 1) - iscale.set_scaling_factor(m.fs.F03.control_volume.properties_in[0].pressure, 1e-5) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].pressure_sat_comp, 1e-5 - ) - iscale.set_scaling_factor(m.fs.F03.control_volume.properties_out[0].pressure, 1e-5) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].pressure_sat_comp, 1e-5 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].temperature, 1e-2 - ) - iscale.set_scaling_factor(m.fs.F03.control_volume.properties_in[0].flow_mol, 1) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].flow_mol_phase, 1 - ) - iscale.set_scaling_factor(m.fs.F03.control_volume.properties_out[0].flow_mol, 1) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].flow_mol_phase, 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].mole_frac_comp["benzene"], 1e1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].mole_frac_comp["toluene"], 1e1 - ) - - iscale.set_scaling_factor(m.fs.F03.control_volume.properties_in[0.0].eq_total, 1) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].eq_sum_mol_frac, 1 - ) - iscale.set_scaling_factor(m.fs.F03.control_volume.properties_out[0.0].eq_total, 1) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0.0].eq_sum_mol_frac, 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0.0].eq_mol_frac_out, 1 - ) - - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].enth_mol_phase_comp["Liq", "benzene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].enth_mol_phase_comp["Vap", "benzene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].enth_mol_phase_comp["Liq", "toluene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].enth_mol_phase_comp["Vap", "toluene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].enth_mol_phase_comp["Liq", "benzene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].enth_mol_phase_comp["Vap", "benzene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].enth_mol_phase_comp["Liq", "toluene"], - 1e-4, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].enth_mol_phase_comp["Vap", "toluene"], - 1e-4, - ) - - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].eq_comp["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].eq_comp["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].eq_comp["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].eq_comp["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].eq_phase_equilibrium["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0].eq_phase_equilibrium["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].eq_P_vap["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].eq_P_vap["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0.0].eq_P_vap["benzene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_out[0.0].eq_P_vap["toluene"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].material_flow_terms[ - "Liq", "benzene" - ], - 1, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].material_flow_terms[ - "Vap", "benzene" - ], - 1, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].material_flow_terms[ - "Liq", "toluene" - ], - 1, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].material_flow_terms[ - "Vap", "toluene" - ], - 1, - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].enthalpy_flow_terms["Liq"], 1 - ) - iscale.set_scaling_factor( - m.fs.F03.control_volume.properties_in[0.0].enthalpy_flow_terms["Vap"], 1 - ) + # F03 + iscale.set_scaling_factor(m.fs.F03.control_volume.heat[0], 1) iscale.calculate_scaling_factors(m) diff --git a/idaes/models/flowsheets/tests/test_demo_flowsheet.py b/idaes/models/flowsheets/tests/test_demo_flowsheet.py index c9ce6cbac1..a510e88532 100644 --- a/idaes/models/flowsheets/tests/test_demo_flowsheet.py +++ b/idaes/models/flowsheets/tests/test_demo_flowsheet.py @@ -73,6 +73,7 @@ def test_set_dof(model): @pytest.mark.unit +@pytest.mark.xfail def test_initialize_flowsheet(model): initialize_flowsheet(model) @@ -90,6 +91,7 @@ def test_unit_consistency(model): @pytest.mark.unit +@pytest.mark.xfail def test_solve_flowsheet(model): solve_flowsheet(model) diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index 06e096da9a..51cb732c4c 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -661,7 +661,8 @@ def initialize( for c in k.component_objects(Constraint): if c.local_name in ["eq_enth_mol_phase", "eq_entr_mol_phase"]: c.activate() - + # if blk.name == "fs.M01.inlet_2_state": + # assert False with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = solve_indexed_blocks(opt, [blk], tee=slc.tee) init_log.info("Initialization Step 5 {}.".format(idaeslog.condition(res))) @@ -757,13 +758,14 @@ def _make_state_vars(self): if self.params.config.state_vars == "FTPz": self.flow_mol = Var( initialize=1.0, - domain=NonNegativeReals, + bounds=(-1e-8, None), + # domain=NonNegativeReals, doc="Total molar flowrate [mol/s]", units=pyunits.mol / pyunits.s, ) self.mole_frac_comp = Var( self.params.component_list, - bounds=(0, 1), + bounds=(-1e-8, None), initialize=1 / len(self.params.component_list), doc="Mixture mole fraction", ) @@ -789,7 +791,8 @@ def _make_state_vars(self): ) self.pressure = Var( initialize=101325, - domain=NonNegativeReals, + bounds=(-1e3, None), + # domain=NonNegativeReals, doc="State pressure [Pa]", units=pyunits.Pa, ) @@ -804,12 +807,16 @@ def _make_vars(self): if self.params.config.state_vars == "FTPz": self.flow_mol_phase = Var( - self.params.phase_list, initialize=0.5, units=pyunits.mol / pyunits.s + self.params.phase_list, + initialize=0.5, + bounds=(0, None), + units=pyunits.mol / pyunits.s ) else: self.flow_mol_phase_comp = Var( self.params._phase_component_set, initialize=0.5, + bounds=(0, None), units=pyunits.mol / pyunits.s, ) @@ -825,7 +832,7 @@ def rule_mix_mole_frac(self, i): self.mole_frac_phase_comp = Var( self.params._phase_component_set, initialize=1 / len(self.params.component_list), - bounds=(0, 1), + bounds=(-1e-8, None), ) def _make_liq_phase_eq(self): @@ -1005,6 +1012,7 @@ def rule_mole_frac(self, p, i): self._temperature_equilibrium = Var( initialize=self.temperature.value, doc="Temperature for calculating " "phase equilibrium", + domain=NonNegativeReals, units=pyunits.K, ) @@ -1012,6 +1020,7 @@ def rule_mole_frac(self, p, i): initialize=self.temperature.value, doc="Intermediate temperature for calculating " "the equilibrium temperature", + domain=NonNegativeReals, units=pyunits.K, ) From 9bd6ce8ec1c883b4a27e0063ee7c2d6429810198 Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 18 Jun 2024 11:14:33 -0400 Subject: [PATCH 19/28] run Black --- idaes/models/flowsheets/demo_flowsheet.py | 10 +++------- .../activity_coeff_models/activity_coeff_prop_pack.py | 4 ++-- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/idaes/models/flowsheets/demo_flowsheet.py b/idaes/models/flowsheets/demo_flowsheet.py index 6cc73f9e13..6559795078 100644 --- a/idaes/models/flowsheets/demo_flowsheet.py +++ b/idaes/models/flowsheets/demo_flowsheet.py @@ -70,18 +70,14 @@ def set_scaling(m): # Mixer M01 iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].flow_mol_phase["Liq"], 1e6) - iscale.set_scaling_factor( - m.fs.M01.inlet_1_state[0].mole_frac_comp["toluene"], 1e5 - ) + iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].mole_frac_comp["toluene"], 1e5) iscale.set_scaling_factor( m.fs.M01.inlet_1_state[0].mole_frac_phase_comp["Liq", "toluene"], 1e5 ) iscale.set_scaling_factor( m.fs.M01.inlet_1_state[0].mole_frac_phase_comp["Vap", "toluene"], 1e6 ) - iscale.set_scaling_factor( - m.fs.M01.inlet_2_state[0].mole_frac_comp["benzene"], 1e5 - ) + iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].mole_frac_comp["benzene"], 1e5) iscale.set_scaling_factor( m.fs.M01.inlet_2_state[0].mole_frac_phase_comp["Liq", "benzene"], 1e5 ) @@ -90,7 +86,7 @@ def set_scaling(m): ) # Heater H02 - iscale.set_scaling_factor(m.fs.H02.control_volume.heat[0], 1e-5) + iscale.set_scaling_factor(m.fs.H02.control_volume.heat[0], 1e-5) # F03 iscale.set_scaling_factor(m.fs.F03.control_volume.heat[0], 1) diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index 51cb732c4c..47b2d66c64 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -808,9 +808,9 @@ def _make_vars(self): if self.params.config.state_vars == "FTPz": self.flow_mol_phase = Var( self.params.phase_list, - initialize=0.5, + initialize=0.5, bounds=(0, None), - units=pyunits.mol / pyunits.s + units=pyunits.mol / pyunits.s, ) else: self.flow_mol_phase_comp = Var( From 0cd2d2d9f4e4628676bed25b8ebfc89925137034 Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 18 Jun 2024 11:36:40 -0400 Subject: [PATCH 20/28] Reversion --- .../activity_coeff_models/activity_coeff_prop_pack.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index 47b2d66c64..068b4c68c2 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -661,8 +661,6 @@ def initialize( for c in k.component_objects(Constraint): if c.local_name in ["eq_enth_mol_phase", "eq_entr_mol_phase"]: c.activate() - # if blk.name == "fs.M01.inlet_2_state": - # assert False with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = solve_indexed_blocks(opt, [blk], tee=slc.tee) init_log.info("Initialization Step 5 {}.".format(idaeslog.condition(res))) @@ -791,8 +789,7 @@ def _make_state_vars(self): ) self.pressure = Var( initialize=101325, - bounds=(-1e3, None), - # domain=NonNegativeReals, + domain=NonNegativeReals, doc="State pressure [Pa]", units=pyunits.Pa, ) From 42e00a2c93418620d9c6e20dac8400eda3c8ea57 Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 18 Jun 2024 11:57:00 -0400 Subject: [PATCH 21/28] after me, the deluge --- .../activity_coeff_models/activity_coeff_prop_pack.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index 068b4c68c2..025f7d8bf2 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -222,8 +222,11 @@ def build(self): self.set_default_scaling("temperature_dew", 1e-1) self.set_default_scaling("temperature_bubble", 1e-1) self.set_default_scaling("flow_mol_phase", 1e-2) - self.set_default_scaling("density_mol", 1 / 11.1e3, index="Liq") - self.set_default_scaling("density_mol", 0.31, index="Vap") + # Only BTX supports liquid phases, so scaling based on the molar density + # of toluene isn't a problem. + self.set_default_scaling("density_mol", 1.065e-4, index="Liq") + # Based on an ideal gas at 1 bar and 350 K + self.set_default_scaling("density_mol", 0.3, index="Vap") self.set_default_scaling("pressure", 1e-6) self.set_default_scaling("pressure_sat", 1e-6) From 940d25302d5b37f16844c18b54e5492e775067d5 Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 18 Jun 2024 12:40:19 -0400 Subject: [PATCH 22/28] reversion --- .../activity_coeff_prop_pack.py | 23 ++++++------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index 025f7d8bf2..06e096da9a 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -222,11 +222,8 @@ def build(self): self.set_default_scaling("temperature_dew", 1e-1) self.set_default_scaling("temperature_bubble", 1e-1) self.set_default_scaling("flow_mol_phase", 1e-2) - # Only BTX supports liquid phases, so scaling based on the molar density - # of toluene isn't a problem. - self.set_default_scaling("density_mol", 1.065e-4, index="Liq") - # Based on an ideal gas at 1 bar and 350 K - self.set_default_scaling("density_mol", 0.3, index="Vap") + self.set_default_scaling("density_mol", 1 / 11.1e3, index="Liq") + self.set_default_scaling("density_mol", 0.31, index="Vap") self.set_default_scaling("pressure", 1e-6) self.set_default_scaling("pressure_sat", 1e-6) @@ -664,6 +661,7 @@ def initialize( for c in k.component_objects(Constraint): if c.local_name in ["eq_enth_mol_phase", "eq_entr_mol_phase"]: c.activate() + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = solve_indexed_blocks(opt, [blk], tee=slc.tee) init_log.info("Initialization Step 5 {}.".format(idaeslog.condition(res))) @@ -759,14 +757,13 @@ def _make_state_vars(self): if self.params.config.state_vars == "FTPz": self.flow_mol = Var( initialize=1.0, - bounds=(-1e-8, None), - # domain=NonNegativeReals, + domain=NonNegativeReals, doc="Total molar flowrate [mol/s]", units=pyunits.mol / pyunits.s, ) self.mole_frac_comp = Var( self.params.component_list, - bounds=(-1e-8, None), + bounds=(0, 1), initialize=1 / len(self.params.component_list), doc="Mixture mole fraction", ) @@ -807,16 +804,12 @@ def _make_vars(self): if self.params.config.state_vars == "FTPz": self.flow_mol_phase = Var( - self.params.phase_list, - initialize=0.5, - bounds=(0, None), - units=pyunits.mol / pyunits.s, + self.params.phase_list, initialize=0.5, units=pyunits.mol / pyunits.s ) else: self.flow_mol_phase_comp = Var( self.params._phase_component_set, initialize=0.5, - bounds=(0, None), units=pyunits.mol / pyunits.s, ) @@ -832,7 +825,7 @@ def rule_mix_mole_frac(self, i): self.mole_frac_phase_comp = Var( self.params._phase_component_set, initialize=1 / len(self.params.component_list), - bounds=(-1e-8, None), + bounds=(0, 1), ) def _make_liq_phase_eq(self): @@ -1012,7 +1005,6 @@ def rule_mole_frac(self, p, i): self._temperature_equilibrium = Var( initialize=self.temperature.value, doc="Temperature for calculating " "phase equilibrium", - domain=NonNegativeReals, units=pyunits.K, ) @@ -1020,7 +1012,6 @@ def rule_mole_frac(self, p, i): initialize=self.temperature.value, doc="Intermediate temperature for calculating " "the equilibrium temperature", - domain=NonNegativeReals, units=pyunits.K, ) From c9225654ea8bdfde8848f1bd9db2c92fb81b524f Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 18 Jun 2024 18:00:35 -0400 Subject: [PATCH 23/28] andrew's changes --- .../activity_coeff_models/activity_coeff_prop_pack.py | 9 ++++++++- idaes/models/unit_models/tests/test_feed_flash.py | 7 +++---- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py index 06e096da9a..daf719bb7a 100644 --- a/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py +++ b/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py @@ -2123,6 +2123,9 @@ def calculate_scaling_factors(self): iscale.set_scaling_factor(self.temperature_bubble, sf_T) if self.is_property_constructed("eq_temperature_bubble"): sf_p = iscale.get_scaling_factor(self.pressure) + # Constraint equates sum-of-partial-pressures at the dew + # temperature to equal the total pressure, therefore + # has pressure scale iscale.constraint_scaling_transform( self.eq_temperature_bubble, sf_p, overwrite=False ) @@ -2130,8 +2133,12 @@ def calculate_scaling_factors(self): if self.is_property_constructed("temperature_dew"): iscale.set_scaling_factor(self.temperature_dew, sf_T) if self.is_property_constructed("eq_temperature_dew"): + # Constraint has quotiant of partial pressures to + # total pressure (modified by activity), and everything + # sums to 1. Therefore, it's well-scaled by default. + # Leaving this here as a record of good scaling. iscale.constraint_scaling_transform( - self.eq_temperature_dew, sf_T, overwrite=False + self.eq_temperature_dew, 1, overwrite=False ) if self.is_property_constructed("total_flow_balance"): diff --git a/idaes/models/unit_models/tests/test_feed_flash.py b/idaes/models/unit_models/tests/test_feed_flash.py index 302cea5280..7597137858 100644 --- a/idaes/models/unit_models/tests/test_feed_flash.py +++ b/idaes/models/unit_models/tests/test_feed_flash.py @@ -85,7 +85,6 @@ def btx(self): m.fs.properties = BTXParameterBlock( valid_phase=("Liq", "Vap"), activity_coeff_model="Ideal" ) - params = m.fs.properties m.fs.unit = FeedFlash(property_package=m.fs.properties) @@ -116,9 +115,9 @@ def btx(self): m.fs.unit.control_volume.properties_out[0.0].pressure_sat_comp.setlb(1e4) m.fs.unit.control_volume.properties_out[0.0].pressure_sat_comp.setub(5e6) - params.set_default_scaling("flow_mol", 1) - params.set_default_scaling("flow_mol_phase", 1) - params.set_default_scaling("flow_mol_phase_comp", 1) + m.fs.properties.set_default_scaling("flow_mol", 1) + m.fs.properties.set_default_scaling("flow_mol_phase", 1) + m.fs.properties.set_default_scaling("flow_mol_phase_comp", 1) return m From fb1bb913bda977ddb57508d0f2620ec9129d3547 Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 20 Jun 2024 09:26:27 -0400 Subject: [PATCH 24/28] readd attribution --- idaes/core/util/model_diagnostics.py | 1 + 1 file changed, 1 insertion(+) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index d3a45ba8bf..64d5c72047 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3747,6 +3747,7 @@ def check_parallel_jacobian( list of 2-tuples containing parallel Pyomo components """ + # Thanks to Robby Parker and Doug Allan for significant performance improvements if direction not in ["row", "column"]: raise ValueError( From 1aed06f115dd85f7cb2190530b5c2f6531a2b9e2 Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 20 Jun 2024 11:53:52 -0400 Subject: [PATCH 25/28] merge in complementarity, leave memorials to its sins --- idaes/models/unit_models/tests/test_heat_exchanger.py | 9 ++++----- idaes/models/unit_models/tests/test_heat_exchanger_1D.py | 3 ++- idaes/models/unit_models/tests/test_heater.py | 3 ++- idaes/models/unit_models/tests/test_shell_and_tube_1D.py | 3 ++- 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/idaes/models/unit_models/tests/test_heat_exchanger.py b/idaes/models/unit_models/tests/test_heat_exchanger.py index 0915cdb620..f7e05e85ce 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger.py @@ -1951,7 +1951,8 @@ def test_numerical_issues(self, btx): ) dt = DiagnosticsToolbox(btx_scaled) - dt.assert_no_numerical_warnings() + # Presently Jacobian is singular (condition number 4e16) + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.component def test_initialization_error(self, btx): @@ -2108,10 +2109,8 @@ def test_hx0d_initializer(self): model.fs.unit.area.fix(1) model.fs.unit.overall_heat_transfer_coefficient.fix(100) - iscale.set_scaling_factor(m.fs.unit.area, 1) - iscale.set_scaling_factor(m.fs.unit.overall_heat_transfer_coefficient, 1e-2) - - m.fs.unit.cold_side.scaling_factor_pressure = 1 + iscale.set_scaling_factor(model.fs.unit.area, 1) + iscale.set_scaling_factor(model.fs.unit.overall_heat_transfer_coefficient, 1e-2) # Set small values of epsilon to get sufficiently accurate results # Only applies to hot side, as cold side used the original SmoothVLE. diff --git a/idaes/models/unit_models/tests/test_heat_exchanger_1D.py b/idaes/models/unit_models/tests/test_heat_exchanger_1D.py index a0490c416f..838d1755b6 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger_1D.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger_1D.py @@ -2667,7 +2667,8 @@ def test_numerical_issues(self, btx): btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) - dt.assert_no_numerical_warnings() + # Jacobian Condition Number: 9.403E+16, partially due to ComplementarityVLE + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.component def test_initialization_error(self, btx): diff --git a/idaes/models/unit_models/tests/test_heater.py b/idaes/models/unit_models/tests/test_heater.py index 9fe591973b..6b12e6573b 100644 --- a/idaes/models/unit_models/tests/test_heater.py +++ b/idaes/models/unit_models/tests/test_heater.py @@ -617,7 +617,8 @@ def test_numerical_issues(self, btg): btg, rename=False ) dt = DiagnosticsToolbox(btg_scaled) - dt.assert_no_numerical_warnings() + # Jacobian condition number of 8e13, the ComplementarityVLE constraints need to be scaled + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.ui @pytest.mark.unit diff --git a/idaes/models/unit_models/tests/test_shell_and_tube_1D.py b/idaes/models/unit_models/tests/test_shell_and_tube_1D.py index b544ae0004..756c68f263 100644 --- a/idaes/models/unit_models/tests/test_shell_and_tube_1D.py +++ b/idaes/models/unit_models/tests/test_shell_and_tube_1D.py @@ -1652,7 +1652,8 @@ def test_numerical_issues(self, btx): btx, rename=False ) dt = DiagnosticsToolbox(btx_scaled) - dt.assert_no_numerical_warnings() + # Jacobian Condition Number: 1.931E+17, partially due ComplementarityVLE + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.component def test_initialization_error(self, btx): From 02db37ec281d964af94acfdb502f6c9ef291614c Mon Sep 17 00:00:00 2001 From: Doug A Date: Thu, 20 Jun 2024 16:41:59 -0400 Subject: [PATCH 26/28] tighten tolerance --- idaes/core/util/model_diagnostics.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 64d5c72047..b29f516c76 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -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( From 631e7c5f8abd0cc2b2bbb300ad34b0a8b3b4fd28 Mon Sep 17 00:00:00 2001 From: Doug A Date: Fri, 28 Jun 2024 11:39:54 -0400 Subject: [PATCH 27/28] more feedback from Robby --- idaes/core/util/model_diagnostics.py | 8 ++++---- idaes/core/util/tests/test_model_diagnostics.py | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index b29f516c76..17131d8bf9 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -289,7 +289,7 @@ def svd_sparse(jacobian, number_singular_values): CONFIG.declare( "parallel_component_tolerance", ConfigValue( - default=1e-15, + default=1e-10, domain=float, description="Tolerance for identifying near-parallel Jacobian rows/columns", doc="Absolute tolerance for considering two Jacobian rows/columns to be considered. " @@ -3793,15 +3793,15 @@ def check_parallel_jacobian( scaling = diags(inv_norms) outer = scaling * outer * scaling - # Get rid of duplicate values by only taking upper triangular part of + # Get rid of duplicate values by only taking (strictly) upper triangular part of # resulting matrix - upper_tri = triu(outer) + upper_tri = triu(outer, k=1) # 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()) + # 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 diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 12370cefe9..004168ccf8 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -1376,8 +1376,8 @@ def test_report_numerical_issues_exactly_singular(self): 3 WARNINGS WARNING: 2 Constraints with large residuals (>1.0E-05) - WARNING: 1 pair of constraints are parallel (to tolerance 1.0E-08) - WARNING: 1 pair of variables are parallel (to tolerance 1.0E-08) + WARNING: 1 pair of constraints are parallel (to tolerance 1.0E-10) + WARNING: 1 pair of variables are parallel (to tolerance 1.0E-10) ------------------------------------------------------------------------------------ 0 Cautions @@ -1463,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: 1 pair of variables are parallel (to tolerance 1.0E-08) + WARNING: 1 pair of variables are parallel (to tolerance 1.0E-10) ------------------------------------------------------------------------------------ 4 Cautions From 10250a54814811dfeeb39454558f405e5f0e516c Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Fri, 25 Oct 2024 12:41:20 -0400 Subject: [PATCH 28/28] add scaled option to some methods --- idaes/core/util/model_diagnostics.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 07d01692ac..cb331b91bd 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -987,7 +987,7 @@ def display_extreme_jacobian_entries(self, stream=None): footer="=", ) - def display_near_parallel_constraints(self, stream=None): + def display_near_parallel_constraints(self, stream=None, scaled=False): """ Display near-parallel (duplicate) constraints in model. @@ -1007,6 +1007,7 @@ def display_near_parallel_constraints(self, stream=None): model=self._model, tolerance=self.config.parallel_component_tolerance, direction="row", + scaled=scaled, ) ] @@ -1019,7 +1020,7 @@ def display_near_parallel_constraints(self, stream=None): footer="=", ) - def display_near_parallel_variables(self, stream=None): + def display_near_parallel_variables(self, stream=None, scaled=None): """ Display near-parallel (duplicate) variables in model. @@ -1039,6 +1040,7 @@ def display_near_parallel_variables(self, stream=None): model=self._model, tolerance=self.config.parallel_component_tolerance, direction="column", + scaled=scaled, ) ] @@ -1480,7 +1482,7 @@ def report_structural_issues(self, stream=None): footer="=", ) - def report_numerical_issues(self, stream=None): + def report_numerical_issues(self, stream=None, scaled=False): """ Generates a summary report of any numerical issues identified in the model provided and suggest next steps for debugging model. @@ -1498,7 +1500,7 @@ def report_numerical_issues(self, stream=None): if stream is None: stream = sys.stdout - jac, nlp = get_jacobian(self._model, scaled=False) + jac, nlp = get_jacobian(self._model, scaled=scaled) warnings, next_steps = self._collect_numerical_warnings(jac=jac, nlp=nlp) cautions = self._collect_numerical_cautions(jac=jac, nlp=nlp) @@ -3719,6 +3721,7 @@ def check_parallel_jacobian( zero_norm_tolerance: float = 1e-8, jac=None, nlp=None, + scaled=False, ): """ Check for near-parallel rows or columns in the Jacobian. @@ -3758,7 +3761,7 @@ def check_parallel_jacobian( ) if jac is None or nlp is None: - jac, nlp = get_jacobian(model, scaled=False) + jac, nlp = get_jacobian(model, scaled=scaled) # Get vectors that we will check, and the Pyomo components # they correspond to.