From 3c323b14ebb78685dc784006a8a117ff1cb0e731 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 19 Dec 2025 23:51:37 +0100 Subject: [PATCH 01/12] regularize matrices to fix one-time orthogonalization error --- Common/src/linear_algebra/CSysMatrix.cpp | 50 +++++++++++++++++++++++- 1 file changed, 48 insertions(+), 2 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index dab97682f499..17b82772800d 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -505,9 +505,23 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v #define A(I, J) matrix[(I)*nVar + (J)] /*--- Transform system in Upper Matrix ---*/ + + /*--- Regularization epsilon to prevent divide-by-zero ---*/ + //constexpr ScalarType eps = 1e-12; + for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { + + /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ + if (std::abs(A(jVar, jVar)) < EPS) { + ScalarType sign = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); + A(jVar, jVar) = sign * EPS; + std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar + << ") to " << A(jVar, jVar) << std::endl; + } + ScalarType weight = A(iVar, jVar) / A(jVar, jVar); + for (auto kVar = jVar; kVar < nVar; kVar++) A(iVar, kVar) -= weight * A(jVar, kVar); vec[iVar] -= weight * vec[jVar]; } @@ -517,7 +531,17 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v for (auto iVar = nVar; iVar > 0ul;) { iVar--; // unsigned type for (auto jVar = iVar + 1; jVar < nVar; jVar++) vec[iVar] -= A(iVar, jVar) * vec[jVar]; + + /*--- Regularize diagonal if too small ---*/ + if (std::abs(A(iVar, iVar)) < EPS) { + ScalarType sign = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); + A(iVar, iVar) = sign * EPS; + std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar + << ") to " << A(iVar, iVar) << std::endl; + } + vec[iVar] /= A(iVar, iVar); + } #undef A #endif @@ -547,11 +571,22 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver #else #define A(I, J) matrix[(I)*nVar + (J)] + /*--- Regularization epsilon to prevent divide-by-zero ---*/ + //constexpr ScalarType eps = 1e-12; + /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { - ScalarType weight = A(iVar, jVar) / A(jVar, jVar); + /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ + if (std::abs(A(jVar, jVar)) < EPS) { + ScalarType sign = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); + A(jVar, jVar) = sign * EPS; + std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar + << ") to " << A(jVar, jVar) << " at iVar=" << iVar << std::endl; + } + + ScalarType weight = A(iVar, jVar) / A(jVar, jVar); for (auto kVar = jVar; kVar < nVar; kVar++) A(iVar, kVar) -= weight * A(jVar, kVar); /*--- at this stage M is lower triangular so not all cols need updating ---*/ @@ -565,7 +600,18 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver for (auto jVar = iVar + 1; jVar < nVar; jVar++) for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) -= A(iVar, jVar) * M(jVar, kVar); - for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) /= A(iVar, iVar); + /*--- Regularize diagonal if too small ---*/ + if (std::abs(A(iVar, iVar)) < EPS) { + ScalarType sign = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); + A(iVar, iVar) = sign * EPS; + std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar + << ") to " << A(iVar, iVar) << std::endl; + } + + for (auto kVar = 0ul; kVar < nVar; kVar++) { + M(iVar, kVar) /= A(iVar, iVar); + + } } #undef A #endif From 10fada816a06da283cc0d7f21eb4c06caa018146 Mon Sep 17 00:00:00 2001 From: Nijso Date: Fri, 19 Dec 2025 23:56:30 +0100 Subject: [PATCH 02/12] Apply suggestion from @bigfooted --- Common/src/linear_algebra/CSysMatrix.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 17b82772800d..ed07dd3a0d2b 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -571,9 +571,6 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver #else #define A(I, J) matrix[(I)*nVar + (J)] - /*--- Regularization epsilon to prevent divide-by-zero ---*/ - //constexpr ScalarType eps = 1e-12; - /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { From ea1e4373e63bebaa69c535e5b8a4554812f24476 Mon Sep 17 00:00:00 2001 From: Nijso Date: Fri, 19 Dec 2025 23:56:49 +0100 Subject: [PATCH 03/12] Apply suggestion from @bigfooted --- Common/src/linear_algebra/CSysMatrix.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index ed07dd3a0d2b..33ffc30a09e3 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -505,10 +505,6 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v #define A(I, J) matrix[(I)*nVar + (J)] /*--- Transform system in Upper Matrix ---*/ - - /*--- Regularization epsilon to prevent divide-by-zero ---*/ - //constexpr ScalarType eps = 1e-12; - for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { From b202ee0f902342628183875abef40c8bc46ad925 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 20 Dec 2025 06:53:01 +0100 Subject: [PATCH 04/12] precommit --- Common/src/linear_algebra/CSysMatrix.cpp | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 17b82772800d..56e0173db43a 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -507,17 +507,16 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v /*--- Transform system in Upper Matrix ---*/ /*--- Regularization epsilon to prevent divide-by-zero ---*/ - //constexpr ScalarType eps = 1e-12; + // constexpr ScalarType eps = 1e-12; for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { - /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < EPS) { ScalarType sign = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); A(jVar, jVar) = sign * EPS; - std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar - << ") to " << A(jVar, jVar) << std::endl; + std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar << ") to " + << A(jVar, jVar) << std::endl; } ScalarType weight = A(iVar, jVar) / A(jVar, jVar); @@ -536,12 +535,11 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v if (std::abs(A(iVar, iVar)) < EPS) { ScalarType sign = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); A(iVar, iVar) = sign * EPS; - std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar - << ") to " << A(iVar, iVar) << std::endl; + std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " + << A(iVar, iVar) << std::endl; } vec[iVar] /= A(iVar, iVar); - } #undef A #endif @@ -572,18 +570,17 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver #define A(I, J) matrix[(I)*nVar + (J)] /*--- Regularization epsilon to prevent divide-by-zero ---*/ - //constexpr ScalarType eps = 1e-12; + // constexpr ScalarType eps = 1e-12; /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { - /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < EPS) { ScalarType sign = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); A(jVar, jVar) = sign * EPS; - std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar - << ") to " << A(jVar, jVar) << " at iVar=" << iVar << std::endl; + std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to " + << A(jVar, jVar) << " at iVar=" << iVar << std::endl; } ScalarType weight = A(iVar, jVar) / A(jVar, jVar); @@ -604,13 +601,12 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver if (std::abs(A(iVar, iVar)) < EPS) { ScalarType sign = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); A(iVar, iVar) = sign * EPS; - std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar - << ") to " << A(iVar, iVar) << std::endl; + std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " + << A(iVar, iVar) << std::endl; } for (auto kVar = 0ul; kVar < nVar; kVar++) { M(iVar, kVar) /= A(iVar, iVar); - } } #undef A From f07d82619e6c4740db63fdc36730e8c713dca86c Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 20 Dec 2025 20:33:34 +0100 Subject: [PATCH 05/12] codi fix for EPS --- Common/src/linear_algebra/CSysMatrix.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 4d9017c69627..224a72076eb1 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -511,7 +511,7 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < EPS) { ScalarType sign = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); - A(jVar, jVar) = sign * EPS; + A(jVar, jVar) = ScalarType(sign * EPS); std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << std::endl; } @@ -531,7 +531,7 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v /*--- Regularize diagonal if too small ---*/ if (std::abs(A(iVar, iVar)) < EPS) { ScalarType sign = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); - A(iVar, iVar) = sign * EPS; + A(iVar, iVar) = ScalarType(sign * EPS); std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " << A(iVar, iVar) << std::endl; } @@ -572,7 +572,7 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < EPS) { ScalarType sign = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); - A(jVar, jVar) = sign * EPS; + A(jVar, jVar) = ScalarType(sign * EPS); std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << " at iVar=" << iVar << std::endl; } @@ -594,7 +594,7 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver /*--- Regularize diagonal if too small ---*/ if (std::abs(A(iVar, iVar)) < EPS) { ScalarType sign = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); - A(iVar, iVar) = sign * EPS; + A(iVar, iVar) = ScalarType(sign * EPS); std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " << A(iVar, iVar) << std::endl; } From f51ccebf1da4ad5bbdc100c0b2883201ae2ac119 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 20 Dec 2025 22:42:44 +0100 Subject: [PATCH 06/12] codi fix for EPS --- Common/src/linear_algebra/CSysMatrix.cpp | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 224a72076eb1..95c0d0f9b99d 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -510,8 +510,7 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < EPS) { - ScalarType sign = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); - A(jVar, jVar) = ScalarType(sign * EPS); + A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(EPS) : ScalarType(-EPS); std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << std::endl; } @@ -530,8 +529,7 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v /*--- Regularize diagonal if too small ---*/ if (std::abs(A(iVar, iVar)) < EPS) { - ScalarType sign = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); - A(iVar, iVar) = ScalarType(sign * EPS); + A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(EPS) : ScalarType(-EPS); std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " << A(iVar, iVar) << std::endl; } @@ -571,8 +569,7 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < EPS) { - ScalarType sign = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); - A(jVar, jVar) = ScalarType(sign * EPS); + A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(EPS) : ScalarType(-EPS); std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << " at iVar=" << iVar << std::endl; } @@ -593,8 +590,7 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver /*--- Regularize diagonal if too small ---*/ if (std::abs(A(iVar, iVar)) < EPS) { - ScalarType sign = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(1) : ScalarType(-1); - A(iVar, iVar) = ScalarType(sign * EPS); + A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(EPS) : ScalarType(-EPS); std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " << A(iVar, iVar) << std::endl; } From 9d375620c9bc7d530e0cc08e3748907694f30409 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 20 Dec 2025 23:28:11 +0100 Subject: [PATCH 07/12] codi fix for EPS --- Common/src/linear_algebra/CSysMatrix.cpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 95c0d0f9b99d..2b1173b4ae2f 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -504,13 +504,16 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v #else #define A(I, J) matrix[(I)*nVar + (J)] + /*--- Regularization epsilon to prevent divide-by-zero ---*/ + constexpr passivedouble eps = 1e-12; + /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ - if (std::abs(A(jVar, jVar)) < EPS) { - A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(EPS) : ScalarType(-EPS); + if (std::abs(A(jVar, jVar)) < eps) { + A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << std::endl; } @@ -528,8 +531,8 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v for (auto jVar = iVar + 1; jVar < nVar; jVar++) vec[iVar] -= A(iVar, jVar) * vec[jVar]; /*--- Regularize diagonal if too small ---*/ - if (std::abs(A(iVar, iVar)) < EPS) { - A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(EPS) : ScalarType(-EPS); + if (std::abs(A(iVar, iVar)) < eps) { + A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " << A(iVar, iVar) << std::endl; } @@ -564,12 +567,15 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver #else #define A(I, J) matrix[(I)*nVar + (J)] + /*--- Regularization epsilon to prevent divide-by-zero ---*/ + constexpr passivedouble eps = 1e-12; + /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ - if (std::abs(A(jVar, jVar)) < EPS) { - A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(EPS) : ScalarType(-EPS); + if (std::abs(A(jVar, jVar)) < eps) { + A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << " at iVar=" << iVar << std::endl; } @@ -589,8 +595,8 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) -= A(iVar, jVar) * M(jVar, kVar); /*--- Regularize diagonal if too small ---*/ - if (std::abs(A(iVar, iVar)) < EPS) { - A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(EPS) : ScalarType(-EPS); + if (std::abs(A(iVar, iVar)) < eps) { + A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " << A(iVar, iVar) << std::endl; } From 04f0681517b3a8cb96a274a6c7ed632b271f44a8 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 21 Dec 2025 20:14:14 +0100 Subject: [PATCH 08/12] eps change --- Common/src/linear_algebra/CSysMatrix.cpp | 15 +- .../turbulent_premixed_psi/QUICK_FIX_GUIDE.md | 124 ++++++++++++ .../turbulent_premixed_psi/SOLUTION_NOTES.md | 128 ++++++++++++ .../turbulent_premixed_psi/VISUALIZATION.md | 187 ++++++++++++++++++ 4 files changed, 448 insertions(+), 6 deletions(-) create mode 100644 TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md create mode 100644 TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md create mode 100644 TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 2b1173b4ae2f..2b20512837a8 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -505,7 +505,8 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v #define A(I, J) matrix[(I)*nVar + (J)] /*--- Regularization epsilon to prevent divide-by-zero ---*/ - constexpr passivedouble eps = 1e-12; + //constexpr passivedouble eps = 1e-12; + const su2double eps = 1e-12; /*--- Transform system in Upper Matrix ---*/ @@ -513,7 +514,8 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < eps) { - A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); + //A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); + A(jVar, jVar) = (A(jVar, jVar) >= 0.0 ? eps : -eps); std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << std::endl; } @@ -532,7 +534,8 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v /*--- Regularize diagonal if too small ---*/ if (std::abs(A(iVar, iVar)) < eps) { - A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); + //A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); + A(iVar, iVar) = (A(iVar, iVar) >= 0.0 ? eps : -eps); std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " << A(iVar, iVar) << std::endl; } @@ -568,14 +571,14 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver #define A(I, J) matrix[(I)*nVar + (J)] /*--- Regularization epsilon to prevent divide-by-zero ---*/ - constexpr passivedouble eps = 1e-12; + const su2double eps = 1e-12; /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < eps) { - A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); + A(jVar, jVar) = (A(jVar, jVar) >= 0.0 ? eps : -eps); std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << " at iVar=" << iVar << std::endl; } @@ -596,7 +599,7 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver /*--- Regularize diagonal if too small ---*/ if (std::abs(A(iVar, iVar)) < eps) { - A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); + A(iVar, iVar) = (A(iVar, iVar) >= 0.0 ? eps : -eps); std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " << A(iVar, iVar) << std::endl; } diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md b/TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md new file mode 100644 index 000000000000..66cbad7f19e9 --- /dev/null +++ b/TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md @@ -0,0 +1,124 @@ +# Quick Fix Summary: Residual Stagnation Solution + +## The Problem +``` +Species equation → converges ✓ +Update T from species → T = (1-c)*Tu + c*Tf +Flow equations → STAGNATE ✗ +``` + +**Why?** Density wasn't updated when temperature changed! + +## The Fix + +### BEFORE (Broken): +```python +def update_temperature(SU2Driver, iPoint): + C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0) + T = Tu*(1-C) + Tf*C + Cp = SU2Driver.Primitives()(iPoint,iCp) + SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iTEMP,Cp*T) + # ⚠️ Density not updated! ρ still has old value +``` + +### AFTER (Fixed): +```python +def update_temperature(SU2Driver, iPoint): + C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0) + T_target = Tu*(1-C) + Tf*C + T_old = SU2Driver.Primitives()(iPoint,iTEMPERATURE) + + # ✓ Under-relaxation for stability + T_new = RELAX_FACTOR * T_target + (1.0 - RELAX_FACTOR) * T_old + + # ✓ Compute gas constant from old state + R_gas = P / (rho_old * T_old) + + # ✓ Update density consistently: ρ = P/(RT) + rho_new = P / (R_gas * T_new) + + # ✓ Update both temperature and density + SU2Driver.Solution(iFLOWSOLVER).Set(iPoint, iTEMP, T_new) + SU2Driver.Primitives().Set(iPoint, iDENSITY, rho_new) +``` + +## Main Loop Sequence + +### BEFORE: +```python +driver.Preprocess() +driver.Run() +# Set source terms (WRONG: after Run!) +for i in nodes: Source.Set(i, zimont(i)) +# Update temperature (no density update) +for i in nodes: update_temperature(i) +driver.Postprocess() +``` + +### AFTER: +```python +driver.Preprocess() +# ✓ Set source terms BEFORE Run +for i in nodes: Source.Set(i, zimont(i)) +# ✓ Run solvers +driver.Run() +# ✓ Update temperature AND density AFTER Run +for i in nodes: update_temperature(i) +# ✓ Monitor coupling +print(f"Max ΔT: {max_delta_T}, Max Δρ/ρ: {max_delta_rho}") +driver.Postprocess() +``` + +## Key Parameters + +```python +# Add at top of file: +RELAX_FACTOR = 0.7 # Adjust 0.5-0.8 for stability vs speed +``` + +## Expected Behavior + +| Before Fix | After Fix | +|------------|-----------| +| Species: ✓ converging | Species: ✓ converging | +| Pressure: ✗ stagnating | Pressure: ✓ converging | +| Velocity: ✗ stagnating | Velocity: ✓ converging | +| ρ inconsistent with T | ρ = P/(RT) ✓ | + +## Physics + +**In combustion:** +- c: 0 → 1 (progress) +- T: 673K → 1777K ↑ +- ρ: 2.52 → ~0.95 ↓ (MUST decrease!) + +**If ρ not updated:** +- Continuity: wrong mass fluxes +- Momentum: wrong inertia +- → Residuals can't converge! + +## Monitoring + +Look for in output: +``` +Source term - Max: 1.23e+02, Min: 0.00e+00, Avg: 3.45e+01 +Max temperature change: 5.23e+01 K +Max relative density change: 2.15e-02 +``` + +## Troubleshooting + +**Still not converging?** +1. Reduce `RELAX_FACTOR` to 0.5 +2. Reduce CFL number in config +3. Check source term isn't too strong +4. Verify mass conservation: Σ(source terms) ≈ 0 + +**Oscillations?** +1. Reduce `RELAX_FACTOR` +2. Check grid quality in flame region +3. Ensure boundary conditions are correct + +## One-Line Summary + +**Always update density when temperature changes: `ρ_new = P/(R*T_new)`** diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md b/TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md new file mode 100644 index 000000000000..bcd29e878fdd --- /dev/null +++ b/TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md @@ -0,0 +1,128 @@ +# Solution for Residual Stagnation in Turbulent Premixed Combustion + +## Problem Description + +When using the Python wrapper to add species transport with strong source terms and algebraically updating temperature/enthalpy from species composition, the pressure and velocity residuals stagnate in regions with high source terms while the species residuals converge properly. + +## Root Cause + +The issue occurs because: + +1. **Species solver** converges and updates the progress variable `c` +2. **Temperature is updated algebraically**: `T = (1-c)*Tu + c*Tf` +3. **Density is NOT updated** to match the new temperature +4. **Thermodynamic inconsistency**: The ideal gas law `ρ = P/(RT)` is violated +5. **Flow equations see inconsistent state**: Continuity and momentum equations operate with outdated density + +This breaks the coupling between species, energy, and momentum equations, causing residual stagnation. + +## Solution Implemented (Solution #1) + +### Key Changes in `run.py`: + +#### 1. **Updated `update_temperature()` function** (Lines ~90-130) + - Now computes new density when temperature changes: `ρ_new = P/(R*T_new)` + - Updates both temperature AND density in the primitive variables + - Added under-relaxation factor for stability: `T_new = ω*T_target + (1-ω)*T_old` + - Maintains thermodynamic consistency with ideal gas law + +#### 2. **Improved iteration sequence** (Lines ~250-290) + - **STEP 1**: Set species source terms + - **STEP 2**: Run all solvers (species equation solved) + - **STEP 3**: Update temperature AND density from new species field + - **STEP 4**: Postprocess and update + - Ensures proper coupling between species and flow equations + +#### 3. **Added under-relaxation** (Line ~58) + - `RELAX_FACTOR = 0.7` prevents oscillations + - Can be tuned between 0.5-0.8 depending on source term strength + - Lower values = more stable but slower convergence + +#### 4. **Added coupling diagnostics** (Lines ~270-280) + - Monitors `max_delta_T` and `max_delta_rho` + - Helps identify coupling issues + - Shows when thermodynamic updates are stabilizing + +#### 5. **Added source term monitoring** (Lines ~185-215, ~260-265) + - Tracks max/min/avg source term values + - Helps diagnose convergence problems + - Important for checking mass conservation + +## Physical Interpretation + +For a premixed flame: +- Progress variable goes from `c=0` (unburnt) to `c=1` (burnt) +- Temperature increases: `T: 673K → 1777K` +- Density MUST decrease (isobaric): `ρ ∝ 1/T` +- Without density update: continuity equation sees wrong mass fluxes +- Result: momentum residuals can't converge + +## Expected Results + +After this fix: +- ✅ Species residuals converge (as before) +- ✅ Temperature correctly updates from species +- ✅ **Density correctly updates from temperature** +- ✅ **Pressure residuals converge** (no longer stagnate) +- ✅ **Velocity residuals converge** (no longer stagnate) +- ✅ Thermodynamic consistency maintained + +## Tuning Parameters + +If convergence is still slow, adjust: + +### 1. **Under-relaxation factor** (Line 58) +```python +RELAX_FACTOR = 0.5-0.8 # Lower = more stable, slower +``` + +### 2. **CFL number in config file** +``` +CFL_NUMBER= 10.0 # Reduce if needed +CFL_ADAPT= YES +CFL_ADAPT_PARAM= ( 0.3, 0.5, 10.0, 100.0, 1e-6 ) +``` + +### 3. **Number of inner iterations** +```python +for inner_iter in range(10): # Increase for more coupling iterations +``` + +## Additional Recommendations + +1. **Monitor residual ratios**: Watch for regions where flow residuals are high relative to species residuals + +2. **Check mass conservation**: The source term should not create/destroy mass: + ``` + ∂(ρc)/∂t + ∇·(ρcv) = S_c + ``` + where `S_c` is per unit volume + +3. **Verify ideal gas law**: At any point, check if `P = ρRT` holds after updates + +4. **Use implicit time stepping**: Helps with stiff source terms (already in config) + +5. **Consider local time stepping**: Can help in regions with varying source terms + +## Technical Details + +### Incompressible Flow Considerations +This case uses `INC.FLOW` (incompressible solver), which: +- Solves for pressure directly (not density-based) +- Assumes low Mach number +- Still requires density updates for variable density flows (like combustion!) +- Energy equation is decoupled but affects momentum through density + +### Variable Density Incompressible Flow +Even though it's "incompressible", combustion creates variable density: +- ∇·(ρv) = 0 (mass conservation, not ∇·v = 0) +- ρ varies due to temperature, not pressure waves +- Density MUST be updated when temperature changes + +## References +- See main SU2 documentation for species transport +- Zimont turbulent flame speed closure: Zimont, V.L., 2000. "Gas premixed combustion at high turbulence." +- Variable density incompressible flow formulation + +## Contact +For issues or questions about this fix, refer to the SU2 forums or GitHub issues. diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md b/TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md new file mode 100644 index 000000000000..f7f1f9c20870 --- /dev/null +++ b/TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md @@ -0,0 +1,187 @@ +# Visualization of the Coupling Problem and Solution + +## The Coupling Problem + +``` +┌─────────────────────────────────────────────────────────────┐ +│ ITERATION N │ +└─────────────────────────────────────────────────────────────┘ + +WITHOUT FIX (BROKEN): +━━━━━━━━━━━━━━━━━━━━━ +┌─────────────┐ +│ Species Eq │ ∂(ρc)/∂t + ∇·(ρcv) = S_source +│ │ → Solves for c with source term +│ c: 0.3→0.5 │ → Converges well ✓ +└─────────────┘ + ↓ +┌─────────────┐ +│ Update T │ T = (1-c)*Tu + c*Tf +│ T: 900→1100 │ → Temperature updates ✓ +└─────────────┘ + ↓ + ✗ ← NO DENSITY UPDATE! + ↓ +┌─────────────┐ +│ Flow Eqs │ ρ still has OLD value (ρ_old = 2.1) +│ (next iter) │ But T is NEW (T_new = 1100 K) +│ │ → Violates ρ = P/(RT) ✗ +│ Residuals: │ → Momentum sees wrong inertia +│ Pressure ✗ │ → Continuity sees wrong mass flux +│ Velocity ✗ │ → STAGNATION! +└─────────────┘ + + +WITH FIX (CORRECT): +━━━━━━━━━━━━━━━━━━━ +┌─────────────┐ +│ Species Eq │ ∂(ρc)/∂t + ∇·(ρcv) = S_source +│ │ → Solves for c with source term +│ c: 0.3→0.5 │ → Converges well ✓ +└─────────────┘ + ↓ +┌─────────────┐ +│ Update T │ T = (1-c)*Tu + c*Tf +│ T: 900→1100 │ → Temperature updates ✓ +│ │ +│ Update ρ │ ρ = P/(R*T_new) +│ ρ: 2.1→1.6 │ → Density updates ✓ +└─────────────┘ + ↓ + ✓ ← CONSISTENT STATE! + ↓ +┌─────────────┐ +│ Flow Eqs │ ρ and T are consistent +│ (next iter) │ ρ*R*T = P ✓ +│ │ → Momentum sees correct inertia +│ Residuals: │ → Continuity sees correct mass flux +│ Pressure ✓ │ → CONVERGES! +│ Velocity ✓ │ +└─────────────┘ +``` + +## Physical Picture + +``` +UNBURNT REGION (c=0): FLAME REGION (0 Date: Sun, 21 Dec 2025 20:15:00 +0100 Subject: [PATCH 09/12] remove junk file --- .../turbulent_premixed_psi/QUICK_FIX_GUIDE.md | 124 ------------ .../turbulent_premixed_psi/SOLUTION_NOTES.md | 128 ------------ .../turbulent_premixed_psi/VISUALIZATION.md | 187 ------------------ 3 files changed, 439 deletions(-) delete mode 100644 TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md delete mode 100644 TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md delete mode 100644 TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md b/TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md deleted file mode 100644 index 66cbad7f19e9..000000000000 --- a/TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md +++ /dev/null @@ -1,124 +0,0 @@ -# Quick Fix Summary: Residual Stagnation Solution - -## The Problem -``` -Species equation → converges ✓ -Update T from species → T = (1-c)*Tu + c*Tf -Flow equations → STAGNATE ✗ -``` - -**Why?** Density wasn't updated when temperature changed! - -## The Fix - -### BEFORE (Broken): -```python -def update_temperature(SU2Driver, iPoint): - C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0) - T = Tu*(1-C) + Tf*C - Cp = SU2Driver.Primitives()(iPoint,iCp) - SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iTEMP,Cp*T) - # ⚠️ Density not updated! ρ still has old value -``` - -### AFTER (Fixed): -```python -def update_temperature(SU2Driver, iPoint): - C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0) - T_target = Tu*(1-C) + Tf*C - T_old = SU2Driver.Primitives()(iPoint,iTEMPERATURE) - - # ✓ Under-relaxation for stability - T_new = RELAX_FACTOR * T_target + (1.0 - RELAX_FACTOR) * T_old - - # ✓ Compute gas constant from old state - R_gas = P / (rho_old * T_old) - - # ✓ Update density consistently: ρ = P/(RT) - rho_new = P / (R_gas * T_new) - - # ✓ Update both temperature and density - SU2Driver.Solution(iFLOWSOLVER).Set(iPoint, iTEMP, T_new) - SU2Driver.Primitives().Set(iPoint, iDENSITY, rho_new) -``` - -## Main Loop Sequence - -### BEFORE: -```python -driver.Preprocess() -driver.Run() -# Set source terms (WRONG: after Run!) -for i in nodes: Source.Set(i, zimont(i)) -# Update temperature (no density update) -for i in nodes: update_temperature(i) -driver.Postprocess() -``` - -### AFTER: -```python -driver.Preprocess() -# ✓ Set source terms BEFORE Run -for i in nodes: Source.Set(i, zimont(i)) -# ✓ Run solvers -driver.Run() -# ✓ Update temperature AND density AFTER Run -for i in nodes: update_temperature(i) -# ✓ Monitor coupling -print(f"Max ΔT: {max_delta_T}, Max Δρ/ρ: {max_delta_rho}") -driver.Postprocess() -``` - -## Key Parameters - -```python -# Add at top of file: -RELAX_FACTOR = 0.7 # Adjust 0.5-0.8 for stability vs speed -``` - -## Expected Behavior - -| Before Fix | After Fix | -|------------|-----------| -| Species: ✓ converging | Species: ✓ converging | -| Pressure: ✗ stagnating | Pressure: ✓ converging | -| Velocity: ✗ stagnating | Velocity: ✓ converging | -| ρ inconsistent with T | ρ = P/(RT) ✓ | - -## Physics - -**In combustion:** -- c: 0 → 1 (progress) -- T: 673K → 1777K ↑ -- ρ: 2.52 → ~0.95 ↓ (MUST decrease!) - -**If ρ not updated:** -- Continuity: wrong mass fluxes -- Momentum: wrong inertia -- → Residuals can't converge! - -## Monitoring - -Look for in output: -``` -Source term - Max: 1.23e+02, Min: 0.00e+00, Avg: 3.45e+01 -Max temperature change: 5.23e+01 K -Max relative density change: 2.15e-02 -``` - -## Troubleshooting - -**Still not converging?** -1. Reduce `RELAX_FACTOR` to 0.5 -2. Reduce CFL number in config -3. Check source term isn't too strong -4. Verify mass conservation: Σ(source terms) ≈ 0 - -**Oscillations?** -1. Reduce `RELAX_FACTOR` -2. Check grid quality in flame region -3. Ensure boundary conditions are correct - -## One-Line Summary - -**Always update density when temperature changes: `ρ_new = P/(R*T_new)`** diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md b/TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md deleted file mode 100644 index bcd29e878fdd..000000000000 --- a/TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md +++ /dev/null @@ -1,128 +0,0 @@ -# Solution for Residual Stagnation in Turbulent Premixed Combustion - -## Problem Description - -When using the Python wrapper to add species transport with strong source terms and algebraically updating temperature/enthalpy from species composition, the pressure and velocity residuals stagnate in regions with high source terms while the species residuals converge properly. - -## Root Cause - -The issue occurs because: - -1. **Species solver** converges and updates the progress variable `c` -2. **Temperature is updated algebraically**: `T = (1-c)*Tu + c*Tf` -3. **Density is NOT updated** to match the new temperature -4. **Thermodynamic inconsistency**: The ideal gas law `ρ = P/(RT)` is violated -5. **Flow equations see inconsistent state**: Continuity and momentum equations operate with outdated density - -This breaks the coupling between species, energy, and momentum equations, causing residual stagnation. - -## Solution Implemented (Solution #1) - -### Key Changes in `run.py`: - -#### 1. **Updated `update_temperature()` function** (Lines ~90-130) - - Now computes new density when temperature changes: `ρ_new = P/(R*T_new)` - - Updates both temperature AND density in the primitive variables - - Added under-relaxation factor for stability: `T_new = ω*T_target + (1-ω)*T_old` - - Maintains thermodynamic consistency with ideal gas law - -#### 2. **Improved iteration sequence** (Lines ~250-290) - - **STEP 1**: Set species source terms - - **STEP 2**: Run all solvers (species equation solved) - - **STEP 3**: Update temperature AND density from new species field - - **STEP 4**: Postprocess and update - - Ensures proper coupling between species and flow equations - -#### 3. **Added under-relaxation** (Line ~58) - - `RELAX_FACTOR = 0.7` prevents oscillations - - Can be tuned between 0.5-0.8 depending on source term strength - - Lower values = more stable but slower convergence - -#### 4. **Added coupling diagnostics** (Lines ~270-280) - - Monitors `max_delta_T` and `max_delta_rho` - - Helps identify coupling issues - - Shows when thermodynamic updates are stabilizing - -#### 5. **Added source term monitoring** (Lines ~185-215, ~260-265) - - Tracks max/min/avg source term values - - Helps diagnose convergence problems - - Important for checking mass conservation - -## Physical Interpretation - -For a premixed flame: -- Progress variable goes from `c=0` (unburnt) to `c=1` (burnt) -- Temperature increases: `T: 673K → 1777K` -- Density MUST decrease (isobaric): `ρ ∝ 1/T` -- Without density update: continuity equation sees wrong mass fluxes -- Result: momentum residuals can't converge - -## Expected Results - -After this fix: -- ✅ Species residuals converge (as before) -- ✅ Temperature correctly updates from species -- ✅ **Density correctly updates from temperature** -- ✅ **Pressure residuals converge** (no longer stagnate) -- ✅ **Velocity residuals converge** (no longer stagnate) -- ✅ Thermodynamic consistency maintained - -## Tuning Parameters - -If convergence is still slow, adjust: - -### 1. **Under-relaxation factor** (Line 58) -```python -RELAX_FACTOR = 0.5-0.8 # Lower = more stable, slower -``` - -### 2. **CFL number in config file** -``` -CFL_NUMBER= 10.0 # Reduce if needed -CFL_ADAPT= YES -CFL_ADAPT_PARAM= ( 0.3, 0.5, 10.0, 100.0, 1e-6 ) -``` - -### 3. **Number of inner iterations** -```python -for inner_iter in range(10): # Increase for more coupling iterations -``` - -## Additional Recommendations - -1. **Monitor residual ratios**: Watch for regions where flow residuals are high relative to species residuals - -2. **Check mass conservation**: The source term should not create/destroy mass: - ``` - ∂(ρc)/∂t + ∇·(ρcv) = S_c - ``` - where `S_c` is per unit volume - -3. **Verify ideal gas law**: At any point, check if `P = ρRT` holds after updates - -4. **Use implicit time stepping**: Helps with stiff source terms (already in config) - -5. **Consider local time stepping**: Can help in regions with varying source terms - -## Technical Details - -### Incompressible Flow Considerations -This case uses `INC.FLOW` (incompressible solver), which: -- Solves for pressure directly (not density-based) -- Assumes low Mach number -- Still requires density updates for variable density flows (like combustion!) -- Energy equation is decoupled but affects momentum through density - -### Variable Density Incompressible Flow -Even though it's "incompressible", combustion creates variable density: -- ∇·(ρv) = 0 (mass conservation, not ∇·v = 0) -- ρ varies due to temperature, not pressure waves -- Density MUST be updated when temperature changes - -## References -- See main SU2 documentation for species transport -- Zimont turbulent flame speed closure: Zimont, V.L., 2000. "Gas premixed combustion at high turbulence." -- Variable density incompressible flow formulation - -## Contact -For issues or questions about this fix, refer to the SU2 forums or GitHub issues. diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md b/TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md deleted file mode 100644 index f7f1f9c20870..000000000000 --- a/TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md +++ /dev/null @@ -1,187 +0,0 @@ -# Visualization of the Coupling Problem and Solution - -## The Coupling Problem - -``` -┌─────────────────────────────────────────────────────────────┐ -│ ITERATION N │ -└─────────────────────────────────────────────────────────────┘ - -WITHOUT FIX (BROKEN): -━━━━━━━━━━━━━━━━━━━━━ -┌─────────────┐ -│ Species Eq │ ∂(ρc)/∂t + ∇·(ρcv) = S_source -│ │ → Solves for c with source term -│ c: 0.3→0.5 │ → Converges well ✓ -└─────────────┘ - ↓ -┌─────────────┐ -│ Update T │ T = (1-c)*Tu + c*Tf -│ T: 900→1100 │ → Temperature updates ✓ -└─────────────┘ - ↓ - ✗ ← NO DENSITY UPDATE! - ↓ -┌─────────────┐ -│ Flow Eqs │ ρ still has OLD value (ρ_old = 2.1) -│ (next iter) │ But T is NEW (T_new = 1100 K) -│ │ → Violates ρ = P/(RT) ✗ -│ Residuals: │ → Momentum sees wrong inertia -│ Pressure ✗ │ → Continuity sees wrong mass flux -│ Velocity ✗ │ → STAGNATION! -└─────────────┘ - - -WITH FIX (CORRECT): -━━━━━━━━━━━━━━━━━━━ -┌─────────────┐ -│ Species Eq │ ∂(ρc)/∂t + ∇·(ρcv) = S_source -│ │ → Solves for c with source term -│ c: 0.3→0.5 │ → Converges well ✓ -└─────────────┘ - ↓ -┌─────────────┐ -│ Update T │ T = (1-c)*Tu + c*Tf -│ T: 900→1100 │ → Temperature updates ✓ -│ │ -│ Update ρ │ ρ = P/(R*T_new) -│ ρ: 2.1→1.6 │ → Density updates ✓ -└─────────────┘ - ↓ - ✓ ← CONSISTENT STATE! - ↓ -┌─────────────┐ -│ Flow Eqs │ ρ and T are consistent -│ (next iter) │ ρ*R*T = P ✓ -│ │ → Momentum sees correct inertia -│ Residuals: │ → Continuity sees correct mass flux -│ Pressure ✓ │ → CONVERGES! -│ Velocity ✓ │ -└─────────────┘ -``` - -## Physical Picture - -``` -UNBURNT REGION (c=0): FLAME REGION (0 Date: Mon, 22 Dec 2025 23:22:31 +0100 Subject: [PATCH 10/12] Update Common/src/linear_algebra/CSysMatrix.cpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- Common/src/linear_algebra/CSysMatrix.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 2b20512837a8..e48ccf970fe6 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -571,7 +571,7 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver #define A(I, J) matrix[(I)*nVar + (J)] /*--- Regularization epsilon to prevent divide-by-zero ---*/ - const su2double eps = 1e-12; + const float eps = 1e-12; /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { From b379bfdd6e92f84eda4b3cdaa49e350b5f115529 Mon Sep 17 00:00:00 2001 From: Nijso Date: Mon, 22 Dec 2025 23:28:45 +0100 Subject: [PATCH 11/12] Update Common/src/linear_algebra/CSysMatrix.cpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- Common/src/linear_algebra/CSysMatrix.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index e48ccf970fe6..8fbafec22700 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -578,9 +578,11 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ if (std::abs(A(jVar, jVar)) < eps) { - A(jVar, jVar) = (A(jVar, jVar) >= 0.0 ? eps : -eps); - std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to " + A(jVar, jVar) = std::copysign(eps, SU2_TYPE::GetValue(A(jVar, jVar))); +#ifndef NDEBUG + std::cout << "MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to " << A(jVar, jVar) << " at iVar=" << iVar << std::endl; +#endif } ScalarType weight = A(iVar, jVar) / A(jVar, jVar); From 34dc00382c55527c0b456c1d4a56f4687e8960b7 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Mon, 22 Dec 2025 23:40:32 +0100 Subject: [PATCH 12/12] helper function --- Common/src/linear_algebra/CSysMatrix.cpp | 47 ++++++++---------------- 1 file changed, 16 insertions(+), 31 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 8fbafec22700..9af4802ee35a 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -494,6 +494,18 @@ void CSysMatrix::SetValDiagonalZero() { END_SU2_OMP_FOR } +/*--- Helper function to regularize small pivots ---*/ +template +inline void RegularizePivot(ScalarType& pivot, unsigned long row, unsigned long col, const char* context) { + const float eps = 1e-12; + if (std::abs(pivot) < eps) { + pivot = std::copysign(eps, SU2_TYPE::GetValue(pivot)); +#ifndef NDEBUG + std::cout << context << ": Regularized small pivot A(" << row << "," << col << ") to " << pivot << std::endl; +#endif + } +} + template void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* vec) const { #ifdef USE_MKL_LAPACK @@ -504,21 +516,12 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v #else #define A(I, J) matrix[(I)*nVar + (J)] - /*--- Regularization epsilon to prevent divide-by-zero ---*/ - //constexpr passivedouble eps = 1e-12; - const su2double eps = 1e-12; - /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ - if (std::abs(A(jVar, jVar)) < eps) { - //A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); - A(jVar, jVar) = (A(jVar, jVar) >= 0.0 ? eps : -eps); - std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar << ") to " - << A(jVar, jVar) << std::endl; - } + RegularizePivot(A(jVar, jVar), jVar, jVar, "DEBUG Gauss_Elimination"); ScalarType weight = A(iVar, jVar) / A(jVar, jVar); @@ -533,12 +536,7 @@ void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* v for (auto jVar = iVar + 1; jVar < nVar; jVar++) vec[iVar] -= A(iVar, jVar) * vec[jVar]; /*--- Regularize diagonal if too small ---*/ - if (std::abs(A(iVar, iVar)) < eps) { - //A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps); - A(iVar, iVar) = (A(iVar, iVar) >= 0.0 ? eps : -eps); - std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " - << A(iVar, iVar) << std::endl; - } + RegularizePivot(A(iVar, iVar), iVar, iVar, "DEBUG Gauss_Elimination backsubst"); vec[iVar] /= A(iVar, iVar); } @@ -570,20 +568,11 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver #else #define A(I, J) matrix[(I)*nVar + (J)] - /*--- Regularization epsilon to prevent divide-by-zero ---*/ - const float eps = 1e-12; - /*--- Transform system in Upper Matrix ---*/ for (auto iVar = 1ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < iVar; jVar++) { /*--- Regularize pivot if too small to prevent divide-by-zero ---*/ - if (std::abs(A(jVar, jVar)) < eps) { - A(jVar, jVar) = std::copysign(eps, SU2_TYPE::GetValue(A(jVar, jVar))); -#ifndef NDEBUG - std::cout << "MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to " - << A(jVar, jVar) << " at iVar=" << iVar << std::endl; -#endif - } + RegularizePivot(A(jVar, jVar), jVar, jVar, "MatrixInverse"); ScalarType weight = A(iVar, jVar) / A(jVar, jVar); for (auto kVar = jVar; kVar < nVar; kVar++) A(iVar, kVar) -= weight * A(jVar, kVar); @@ -600,11 +589,7 @@ void CSysMatrix::MatrixInverse(ScalarType* matrix, ScalarType* inver for (auto kVar = 0ul; kVar < nVar; kVar++) M(iVar, kVar) -= A(iVar, jVar) * M(jVar, kVar); /*--- Regularize diagonal if too small ---*/ - if (std::abs(A(iVar, iVar)) < eps) { - A(iVar, iVar) = (A(iVar, iVar) >= 0.0 ? eps : -eps); - std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to " - << A(iVar, iVar) << std::endl; - } + RegularizePivot(A(iVar, iVar), iVar, iVar, "DEBUG MatrixInverse backsubst"); for (auto kVar = 0ul; kVar < nVar; kVar++) { M(iVar, kVar) /= A(iVar, iVar);