From 9a2d329e1b2ab456a061a15f49acf62155d46e46 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 29 Jul 2026 19:36:28 +0200 Subject: [PATCH 1/7] fix order --- Common/include/option_structure.hpp | 1 + Common/src/CConfig.cpp | 2 + Common/src/geometry/CMultiGridGeometry.cpp | 230 ++++++++++++--------- config_template.cfg | 5 + 4 files changed, 141 insertions(+), 97 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1ba1bab61fe..5eda82d9fde 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1126,6 +1126,7 @@ struct CMGOptions { su2double MG_Smooth_StagnationTol{0.0}; /*!< \brief Stagnation early exit: stop if current_rms >= prev_rms * tol. 0 = disabled. */ bool MG_Implicit_Lines{false}; /*!< \brief Enable implicit-lines agglomeration from walls. */ unsigned long MG_Implicit_Lines_MaxLength{20}; /*!< \brief Maximum nodes on a wall-normal implicit line (including wall seed). */ + bool MG_Implicit_Lines_Isotropic{false}; /*!< \brief Use isotropic (vs anisotropic) agglomeration along implicit lines. */ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 2c0dce07e68..09d82d8574a 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2069,6 +2069,8 @@ void CConfig::SetConfig_Options() { addBoolOption("MG_IMPLICIT_LINES", MGOptions.MG_Implicit_Lines, false); /*!\brief MG_IMPLICIT_LINES_MAX_LENGTH\n DESCRIPTION: Maximum number of nodes on a wall-normal implicit agglomeration line (including the wall seed node). DEFAULT: 20 \ingroup Config*/ addUnsignedLongOption("MG_IMPLICIT_LINES_MAX_LENGTH", MGOptions.MG_Implicit_Lines_MaxLength, 20); + /*!\brief MG_IMPLICIT_LINES_ISOTROPIC\n DESCRIPTION: Use isotropic agglomeration along implicit lines (4 cells per coarse CV) instead of anisotropic (2 cells per coarse CV). DEFAULT: NO \ingroup Config*/ + addBoolOption("MG_IMPLICIT_LINES_ISOTROPIC", MGOptions.MG_Implicit_Lines_Isotropic, false); /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i). If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/ addDoubleListOption("MG_CFL_SCALING", nMG_CflScaling_p, MG_CflScaling_p); diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 684da742b13..79cd4a6b609 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1286,6 +1286,7 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, const su2double ANGLE_THRESHOLD_DEG = 20.0; /*!< Stop line if direction deviates more than this. */ const unsigned long MAX_LINE_LENGTH = config->GetMGOptions().MG_Implicit_Lines_MaxLength; const su2double cos_threshold = cos(ANGLE_THRESHOLD_DEG * PI_NUMBER / 180.0); + const bool ISOTROPIC = config->GetMGOptions().MG_Implicit_Lines_Isotropic; const unsigned long nPointFine = fine_grid->GetnPoint(); @@ -1381,14 +1382,20 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, if (rank == MASTER_NODE) { cout << "Implicit line agglomeration: detected " << lines.size() << " lines." << endl; + cout << " Mode: " << (ISOTROPIC ? "ISOTROPIC" : "ANISOTROPIC") << endl; } - /*--- Advancing-front greedy pairing with cross-line merging. - * For each pair stage k, process interior positions (1+2k, 1+2k+1). - * When two lines share the same wall-node parent CV, merge their pairs - * into a single 4-child coarse CV. Otherwise create 2-child coarse CVs. ---*/ + /*--- Agglomeration strategy: + * ANISOTROPIC (default): Pair nodes at the SAME distance from wall on DIFFERENT lines. + * Each coarse CV has 2 fine children (from adjacent lines). + * Reduces mesh by factor ~2 normal to wall, preserves resolution along wall. + * + * ISOTROPIC: Group 4 nodes (2 positions × 2 lines) into one coarse CV. + * Each coarse CV has 4 fine children. + * Reduces mesh uniformly by factor ~4 in all directions. + ---*/ vector reserved(nPointFine, 0); - unsigned pair_idx = 0; + unsigned position_idx = 0; while (true) { bool any_work = false; @@ -1399,111 +1406,140 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, for (unsigned long li = 0; li < lines.size(); ++li) { const auto& L = lines[li]; if (L.empty()) continue; - const auto idx2 = 1 + 2 * pair_idx + 1; - if (L.size() <= idx2) continue; // no pair at this stage + if (ISOTROPIC) { + const auto idx2 = 1 + 2 * position_idx + 1; + if (L.size() <= idx2) continue; // no pair at this stage + } else { + if (L.size() <= 1 + position_idx) continue; // no position at this index + } const auto pW = fine_grid->nodes->GetParent_CV(L[0]); parent_to_lines[pW].push_back(li); } vector line_processed(lines.size(), 0); - /*--- A) Cross-line merges: parents with multiple lines ---*/ - for (auto& [parent, line_ids] : parent_to_lines) { - if (line_ids.size() < 2) continue; - - for (size_t k = 0; k + 1 < line_ids.size(); k += 2) { - const auto li1 = line_ids[k]; - const auto li2 = line_ids[k + 1]; - if (line_processed[li1] || line_processed[li2]) continue; - - const auto& L1 = lines[li1]; - const auto& L2 = lines[li2]; - const auto idx1 = 1 + 2 * pair_idx; - const auto idx2 = idx1 + 1; - if (L1.size() <= idx2 || L2.size() <= idx2) continue; - - const auto a = L1[idx1], b = L1[idx2]; - const auto c = L2[idx1], d = L2[idx2]; - - /*--- Skip if any node is already claimed ---*/ - if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b) || - fine_grid->nodes->GetAgglomerate(c) || fine_grid->nodes->GetAgglomerate(d)) - continue; - if (reserved[a] || reserved[b] || reserved[c] || reserved[d]) continue; - - /*--- Geometrical quality check ---*/ - if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config) || - !GeometricalCheck(c, fine_grid, config) || !GeometricalCheck(d, fine_grid, config)) - continue; - - /*--- Guard against duplicate indices ---*/ - if (a == b || a == c || a == d || b == c || b == d || c == d) { - for (auto other_li : line_ids) line_processed[other_li] = 1; - continue; - } - - /*--- Create 4-child coarse CV ---*/ - fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 0, a); - fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 1, b); - fine_grid->nodes->SetParent_CV(c, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 2, c); - fine_grid->nodes->SetParent_CV(d, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 3, d); - nodes->SetnChildren_CV(Index_CoarseCV, 4); - - reserved[a] = reserved[b] = reserved[c] = reserved[d] = 1; - MGQueue_InnerCV.RemoveCV(a); - MGQueue_InnerCV.RemoveCV(b); - MGQueue_InnerCV.RemoveCV(c); - MGQueue_InnerCV.RemoveCV(d); + if (ISOTROPIC) { + /*--- ISOTROPIC MODE: Group 4 children per coarse CV (2 positions × 2 lines) ---*/ + for (auto& [parent, line_ids] : parent_to_lines) { + if (line_ids.size() < 2) continue; + + for (size_t k = 0; k + 1 < line_ids.size(); k += 2) { + const auto li1 = line_ids[k]; + const auto li2 = line_ids[k + 1]; + if (line_processed[li1] || line_processed[li2]) continue; + + const auto& L1 = lines[li1]; + const auto& L2 = lines[li2]; + const auto idx1 = 1 + 2 * position_idx; + const auto idx2 = idx1 + 1; + if (L1.size() <= idx2 || L2.size() <= idx2) continue; + + const auto a = L1[idx1], b = L1[idx2]; + const auto c = L2[idx1], d = L2[idx2]; + + /*--- Skip if any node is already claimed ---*/ + if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b) || + fine_grid->nodes->GetAgglomerate(c) || fine_grid->nodes->GetAgglomerate(d)) + continue; + if (reserved[a] || reserved[b] || reserved[c] || reserved[d]) continue; + + /*--- Geometrical quality check ---*/ + if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config) || + !GeometricalCheck(c, fine_grid, config) || !GeometricalCheck(d, fine_grid, config)) + continue; + + /*--- Guard against duplicate indices ---*/ + if (a == b || a == c || a == d || b == c || b == d || c == d) { + for (auto other_li : line_ids) line_processed[other_li] = 1; + continue; + } - Index_CoarseCV++; - line_processed[li1] = line_processed[li2] = 1; - for (auto other_li : line_ids) - if (other_li != li1 && other_li != li2) line_processed[other_li] = 1; - any_work = true; + /*--- Create 4-child coarse CV (isotropic agglomeration) ---*/ + fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 0, a); + fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 1, b); + fine_grid->nodes->SetParent_CV(c, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 2, c); + fine_grid->nodes->SetParent_CV(d, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 3, d); + nodes->SetnChildren_CV(Index_CoarseCV, 4); + + reserved[a] = reserved[b] = reserved[c] = reserved[d] = 1; + MGQueue_InnerCV.RemoveCV(a); + MGQueue_InnerCV.RemoveCV(b); + MGQueue_InnerCV.RemoveCV(c); + MGQueue_InnerCV.RemoveCV(d); + + Index_CoarseCV++; + line_processed[li1] = line_processed[li2] = 1; + for (auto other_li : line_ids) + if (other_li != li1 && other_li != li2) line_processed[other_li] = 1; + any_work = true; + } + } + } else { + /*--- ANISOTROPIC MODE: Pair nodes at SAME position on DIFFERENT lines ---*/ + for (auto& [parent, line_ids] : parent_to_lines) { + if (line_ids.size() < 2) continue; + + /*--- Pair consecutive lines at the same position ---*/ + for (size_t k = 0; k + 1 < line_ids.size(); k += 2) { + const auto li1 = line_ids[k]; + const auto li2 = line_ids[k + 1]; + if (line_processed[li1] || line_processed[li2]) continue; + + const auto& L1 = lines[li1]; + const auto& L2 = lines[li2]; + const auto pos = 1 + position_idx; + if (L1.size() <= pos || L2.size() <= pos) continue; + + const auto a = L1[pos]; + const auto b = L2[pos]; + + /*--- Skip if any node is already claimed ---*/ + if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b)) continue; + if (reserved[a] || reserved[b]) continue; + + /*--- Geometrical quality check ---*/ + if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config)) continue; + + /*--- Create 2-child coarse CV (anisotropic: same position, different lines) ---*/ + fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 0, a); + fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 1, b); + nodes->SetnChildren_CV(Index_CoarseCV, 2); + + reserved[a] = reserved[b] = 1; + MGQueue_InnerCV.RemoveCV(a); + MGQueue_InnerCV.RemoveCV(b); + + Index_CoarseCV++; + line_processed[li1] = line_processed[li2] = 1; + any_work = true; + } } } - /*--- B) Single-line 2-child merges for remaining lines ---*/ - for (unsigned long li = 0; li < lines.size(); ++li) { - if (line_processed[li]) continue; - const auto& L = lines[li]; - const auto idx1 = 1 + 2 * pair_idx; - const auto idx2 = idx1 + 1; - if (L.size() <= idx2) continue; - - const auto a = L[idx1], b = L[idx2]; - if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b)) continue; - if (reserved[a] || reserved[b]) continue; - if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config)) continue; - - /*--- Create 2-child coarse CV ---*/ - fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 0, a); - fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 1, b); - nodes->SetnChildren_CV(Index_CoarseCV, 2); - - reserved[a] = reserved[b] = 1; - MGQueue_InnerCV.RemoveCV(a); - MGQueue_InnerCV.RemoveCV(b); - - Index_CoarseCV++; - any_work = true; - } - - pair_idx++; + position_idx++; if (!any_work) break; - /*--- Check if any line still has pairs at the next stage ---*/ + /*--- Check if any line still has positions available ---*/ bool any_more = false; - for (const auto& L : lines) { - if (L.size() > 1 + 2 * pair_idx + 1) { - any_more = true; - break; + if (ISOTROPIC) { + for (const auto& L : lines) { + if (L.size() > 1 + 2 * position_idx + 1) { + any_more = true; + break; + } + } + } else { + for (const auto& L : lines) { + if (L.size() > 1 + position_idx) { + any_more = true; + break; + } } } if (!any_more) break; diff --git a/config_template.cfg b/config_template.cfg index a7d357240e0..9e2fe09703d 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1719,6 +1719,11 @@ MG_IMPLICIT_LINES= NO % Maximum nodes on a wall-normal implicit agglomeration line, including the wall seed. % Increase to extend the line deeper into the boundary layer (default 20). MG_IMPLICIT_LINES_MAX_LENGTH= 20 +% +% Use isotropic (vs anisotropic) agglomeration for implicit lines (NO, YES) +% Anisotropic (NO): Pair cells normal to wall (2 cells per coarse CV, reduces mesh ~2x) +% Isotropic (YES): Pair cells in all directions (4 cells per coarse CV, reduces mesh ~4x) +MG_IMPLICIT_LINES_ISOTROPIC= NO % -------------------------- MESH SMOOTHING -----------------------------% % From 06790fb70c60c845d9f0c3deeb2faccd78eb4300 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Thu, 30 Jul 2026 20:18:30 +0200 Subject: [PATCH 2/7] minor implicit line changes --- Common/include/option_structure.hpp | 1 + Common/src/CConfig.cpp | 2 + Common/src/geometry/CMultiGridGeometry.cpp | 464 ++++++++++++++---- .../src/integration/CMultiGridIntegration.cpp | 11 +- config_template.cfg | 4 + 5 files changed, 371 insertions(+), 111 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 5eda82d9fde..45b2a949b0b 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1127,6 +1127,7 @@ struct CMGOptions { bool MG_Implicit_Lines{false}; /*!< \brief Enable implicit-lines agglomeration from walls. */ unsigned long MG_Implicit_Lines_MaxLength{20}; /*!< \brief Maximum nodes on a wall-normal implicit line (including wall seed). */ bool MG_Implicit_Lines_Isotropic{false}; /*!< \brief Use isotropic (vs anisotropic) agglomeration along implicit lines. */ + unsigned long MG_Startup_Iter{100}; /*!< \brief Number of iterations on coarsest mesh during FMG startup phase. */ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 09d82d8574a..a0bcc8f5621 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2071,6 +2071,8 @@ void CConfig::SetConfig_Options() { addUnsignedLongOption("MG_IMPLICIT_LINES_MAX_LENGTH", MGOptions.MG_Implicit_Lines_MaxLength, 20); /*!\brief MG_IMPLICIT_LINES_ISOTROPIC\n DESCRIPTION: Use isotropic agglomeration along implicit lines (4 cells per coarse CV) instead of anisotropic (2 cells per coarse CV). DEFAULT: NO \ingroup Config*/ addBoolOption("MG_IMPLICIT_LINES_ISOTROPIC", MGOptions.MG_Implicit_Lines_Isotropic, false); + /*!\brief MG_STARTUP_ITER\n DESCRIPTION: Number of iterations on the coarsest mesh during Full Multigrid (FMG) startup phase before advancing to finer meshes. DEFAULT: 100 \ingroup Config*/ + addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100); /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i). If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/ addDoubleListOption("MG_CFL_SCALING", nMG_CflScaling_p, MG_CflScaling_p); diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 79cd4a6b609..c9192e17afd 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -312,8 +312,11 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } /*--- Agglomerate high-aspect-ratio interior nodes along implicit lines from walls. ---*/ + unsigned long Index_CoarseCV_before_implicit_lines = Index_CoarseCV; + unsigned long Index_CoarseCV_after_implicit_lines = Index_CoarseCV; if (config->GetMGOptions().MG_Implicit_Lines) { AgglomerateImplicitLines(Index_CoarseCV, fine_grid, config, MGQueue_InnerCV); + Index_CoarseCV_after_implicit_lines = Index_CoarseCV; } /*--- STEP 2: Agglomerate the domain points. ---*/ @@ -428,6 +431,42 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un nPointDomain = Index_CoarseCV; nPoint = nPointDomain; + /*--- DIAGNOSTIC: Check CV child counts after domain agglomeration ---*/ + if (config->GetMGOptions().MG_Implicit_Lines && (rank == MASTER_NODE)) { + unsigned long nCVs_1child = 0, nCVs_2child = 0, nCVs_3child = 0, nCVs_4child = 0, nCVs_other = 0; + unsigned long n_corrupted_implicit_CVs = 0; + for (auto iCV = Index_CoarseCV_before_implicit_lines; iCV < Index_CoarseCV_after_implicit_lines; iCV++) { + const auto nChildren = nodes->GetnChildren_CV(iCV); + if (nChildren == 1) nCVs_1child++; + else if (nChildren == 2) nCVs_2child++; + else if (nChildren == 3) nCVs_3child++; + else if (nChildren == 4) nCVs_4child++; + else nCVs_other++; + + if (nChildren != 2 && !config->GetMGOptions().MG_Implicit_Lines_Isotropic) { + n_corrupted_implicit_CVs++; + if (n_corrupted_implicit_CVs <= 5) { + cout << " CORRUPTION DETECTED in CV " << iCV << ": has " << nChildren << " children (expected 2)" << endl; + cout << " Children nodes: "; + for (unsigned short iChild = 0; iChild < nChildren; iChild++) { + cout << nodes->GetChildren_CV(iCV, iChild); + if (iChild < nChildren - 1) cout << ", "; + } + cout << endl; + } + } + } + if (n_corrupted_implicit_CVs > 0) { + cout << " AFTER DOMAIN AGGLOMERATION: " << n_corrupted_implicit_CVs + << " implicit line CVs were corrupted (child count != 2)" << endl; + cout << " Distribution in implicit line CVs: 1-child=" << nCVs_1child + << ", 2-child=" << nCVs_2child << ", 3-child=" << nCVs_3child + << ", 4-child=" << nCVs_4child; + if (nCVs_other > 0) cout << ", other=" << nCVs_other; + cout << endl; + } + } + /*--- Check that there are no hanging nodes. Detect isolated points (only 1 neighbor), and merge their children CV's with the neighbor. ---*/ @@ -501,6 +540,58 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } } + /*--- Diagnostic: Check if implicit line CVs were corrupted by hanging node correction ---*/ + if (config->GetMGOptions().MG_Implicit_Lines && (rank == MASTER_NODE)) { + unsigned long nCVs_1child = 0, nCVs_2child = 0, nCVs_3child = 0, nCVs_4child = 0, nCVs_other = 0; + unsigned long n_corrupted_after_hanging = 0; + for (auto iCV = Index_CoarseCV_before_implicit_lines; iCV < Index_CoarseCV_after_implicit_lines; iCV++) { + const auto nChildren = nodes->GetnChildren_CV(iCV); + if (nChildren == 1) nCVs_1child++; + else if (nChildren == 2) nCVs_2child++; + else if (nChildren == 3) nCVs_3child++; + else if (nChildren == 4) nCVs_4child++; + else nCVs_other++; + + if (nChildren != 2 && !config->GetMGOptions().MG_Implicit_Lines_Isotropic) { + n_corrupted_after_hanging++; + } + } + if (n_corrupted_after_hanging > 0) { + cout << " AFTER HANGING NODE CORRECTION: " << n_corrupted_after_hanging + << " implicit line CVs corrupted (child count != 2)" << endl; + cout << " Distribution in implicit line CVs: 1-child=" << nCVs_1child + << ", 2-child=" << nCVs_2child << ", 3-child=" << nCVs_3child + << ", 4-child=" << nCVs_4child; + if (nCVs_other > 0) cout << ", other=" << nCVs_other; + cout << endl; + } + } + + /*--- Final summary of all CVs ---*/ + if (config->GetMGOptions().MG_Implicit_Lines && (rank == MASTER_NODE)) { + cout << " Expected ratio: ~2 nodes per CV (actual: " << fixed << setprecision(2) + << (double)fine_grid->GetnPoint() / (double)nPointDomain << ")" << endl; + + unsigned long nCVs_1child = 0, nCVs_2child = 0, nCVs_3child = 0, nCVs_4child = 0, nCVs_other = 0; + for (auto iCV = 0ul; iCV < nPointDomain; iCV++) { + const auto nChildren = nodes->GetnChildren_CV(iCV); + if (nChildren == 1) nCVs_1child++; + else if (nChildren == 2) nCVs_2child++; + else if (nChildren == 3) nCVs_3child++; + else if (nChildren == 4) nCVs_4child++; + else nCVs_other++; + } + cout << " CV distribution: 1-child=" << nCVs_1child << ", 2-child=" << nCVs_2child + << ", 3-child=" << nCVs_3child << ", 4-child=" << nCVs_4child; + if (nCVs_other > 0) cout << ", other=" << nCVs_other; + cout << endl; + + if (nCVs_3child > 0 || (!config->GetMGOptions().MG_Implicit_Lines_Isotropic && nCVs_4child > 0)) { + cout << " WARNING: Detected unexpected CV child counts (3-child=" << nCVs_3child + << ", 4-child=" << nCVs_4child << " in ANISO mode)" << endl; + } + } + /*--- Reset the neighbor information. ---*/ nodes->ResetPoints(); @@ -659,11 +750,9 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un SetGlobal_nPointDomain(Global_nPointCoarse); if (iMesh != MESH_0) { - /*--- Note: CFL at the coarse levels have a large impact on convergence, - this should be rewritten to use adaptive CFL. ---*/ - const su2double Coeff = 1.5; - const su2double CFL = config->GetCFL(iMesh - 1) / Coeff; - config->SetCFL(iMesh, CFL); + /*--- Initialize coarse-level CFL from config. MG_CFL_SCALING will + apply per-level reductions during the multigrid cycle. ---*/ + config->SetCFL(iMesh, config->GetCFL(MESH_0)); } const su2double ratio = su2double(Global_nPointFine) / su2double(Global_nPointCoarse); @@ -1289,6 +1378,9 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, const bool ISOTROPIC = config->GetMGOptions().MG_Implicit_Lines_Isotropic; const unsigned long nPointFine = fine_grid->GetnPoint(); + const unsigned long starting_Index_CoarseCV = Index_CoarseCV; /*--- Track how many CVs we create ---*/ + const bool DEBUG_OUTPUT = (rank == MASTER_NODE); /*--- Enable detailed diagnostic output ---*/ + const unsigned long DEBUG_CV_LIMIT = 20; /*--- Show details for first N CVs ---*/ /*--- Collect implicit lines starting at viscous (no-slip) wall vertices only. * Seeding from non-wall boundaries (farfield, inlet, outlet, symmetry) would @@ -1383,6 +1475,28 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, if (rank == MASTER_NODE) { cout << "Implicit line agglomeration: detected " << lines.size() << " lines." << endl; cout << " Mode: " << (ISOTROPIC ? "ISOTROPIC" : "ANISOTROPIC") << endl; + /*--- Show line length distribution ---*/ + size_t min_len = ULONG_MAX, max_len = 0; + su2double avg_len = 0.0; + for (const auto& L : lines) { + min_len = min(min_len, L.size()); + max_len = max(max_len, L.size()); + avg_len += L.size(); + } + if (!lines.empty()) avg_len /= lines.size(); + cout << " Line lengths: min=" << min_len << ", max=" << max_len << ", avg=" << std::setprecision(1) << std::fixed << avg_len << endl; + + /*--- Show first few lines for debugging ---*/ + cout << " First 5 lines (showing first 4 nodes):" << endl; + for (size_t i = 0; i < min(size_t(5), lines.size()); ++i) { + cout << " Line " << i << " (len=" << lines[i].size() << "): ["; + for (size_t j = 0; j < min(size_t(4), lines[i].size()); ++j) { + if (j > 0) cout << ", "; + cout << lines[i][j]; + } + if (lines[i].size() > 4) cout << ", ..."; + cout << "]" << endl; + } } /*--- Agglomeration strategy: @@ -1399,10 +1513,11 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, while (true) { bool any_work = false; + vector line_processed(lines.size(), 0); - /*--- Build map: wall parent CV -> list of line indices ---*/ - unordered_map> parent_to_lines; - parent_to_lines.reserve(lines.size()); + /*--- Build list of active lines (have nodes at current position) ---*/ + vector active_lines; + active_lines.reserve(lines.size()); for (unsigned long li = 0; li < lines.size(); ++li) { const auto& L = lines[li]; if (L.empty()) continue; @@ -1412,113 +1527,192 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, } else { if (L.size() <= 1 + position_idx) continue; // no position at this index } - const auto pW = fine_grid->nodes->GetParent_CV(L[0]); - parent_to_lines[pW].push_back(li); + active_lines.push_back(li); } - vector line_processed(lines.size(), 0); - if (ISOTROPIC) { - /*--- ISOTROPIC MODE: Group 4 children per coarse CV (2 positions × 2 lines) ---*/ - for (auto& [parent, line_ids] : parent_to_lines) { - if (line_ids.size() < 2) continue; - - for (size_t k = 0; k + 1 < line_ids.size(); k += 2) { - const auto li1 = line_ids[k]; - const auto li2 = line_ids[k + 1]; - if (line_processed[li1] || line_processed[li2]) continue; - - const auto& L1 = lines[li1]; - const auto& L2 = lines[li2]; - const auto idx1 = 1 + 2 * position_idx; - const auto idx2 = idx1 + 1; - if (L1.size() <= idx2 || L2.size() <= idx2) continue; - - const auto a = L1[idx1], b = L1[idx2]; - const auto c = L2[idx1], d = L2[idx2]; - - /*--- Skip if any node is already claimed ---*/ - if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b) || - fine_grid->nodes->GetAgglomerate(c) || fine_grid->nodes->GetAgglomerate(d)) - continue; - if (reserved[a] || reserved[b] || reserved[c] || reserved[d]) continue; - - /*--- Geometrical quality check ---*/ - if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config) || - !GeometricalCheck(c, fine_grid, config) || !GeometricalCheck(d, fine_grid, config)) - continue; - - /*--- Guard against duplicate indices ---*/ - if (a == b || a == c || a == d || b == c || b == d || c == d) { - for (auto other_li : line_ids) line_processed[other_li] = 1; - continue; + /*--- ISOTROPIC MODE: Group 4 children per coarse CV (2 positions × 2 lines) + Use spatial neighbor search to pair adjacent lines. ---*/ + for (auto li1 : active_lines) { + if (line_processed[li1]) continue; + + const auto& L1 = lines[li1]; + const auto idx1 = 1 + 2 * position_idx; + const auto idx2 = idx1 + 1; + if (L1.size() <= idx2) continue; + + const auto a = L1[idx1], b = L1[idx2]; + if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b)) continue; + if (reserved[a] || reserved[b]) continue; + + /*--- Find nearest neighbor line by checking mesh neighbors of node 'a' ---*/ + unsigned long li2_best = std::numeric_limits::max(); + for (auto neighbor_point : fine_grid->nodes->GetPoints(a)) { + /*--- Check if this neighbor belongs to another unprocessed line at same position ---*/ + for (auto li2 : active_lines) { + if (li2 == li1 || line_processed[li2]) continue; + const auto& L2 = lines[li2]; + if (L2.size() <= idx2) continue; + const auto c = L2[idx1]; + if (c == neighbor_point) { + li2_best = li2; + break; + } } + if (li2_best != std::numeric_limits::max()) break; + } + + if (li2_best == std::numeric_limits::max()) continue; - /*--- Create 4-child coarse CV (isotropic agglomeration) ---*/ - fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 0, a); - fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 1, b); - fine_grid->nodes->SetParent_CV(c, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 2, c); - fine_grid->nodes->SetParent_CV(d, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 3, d); - nodes->SetnChildren_CV(Index_CoarseCV, 4); - - reserved[a] = reserved[b] = reserved[c] = reserved[d] = 1; - MGQueue_InnerCV.RemoveCV(a); - MGQueue_InnerCV.RemoveCV(b); - MGQueue_InnerCV.RemoveCV(c); - MGQueue_InnerCV.RemoveCV(d); - - Index_CoarseCV++; - line_processed[li1] = line_processed[li2] = 1; - for (auto other_li : line_ids) - if (other_li != li1 && other_li != li2) line_processed[other_li] = 1; - any_work = true; + const auto& L2 = lines[li2_best]; + const auto c = L2[idx1], d = L2[idx2]; + + /*--- Skip if any node is already claimed ---*/ + if (fine_grid->nodes->GetAgglomerate(c) || fine_grid->nodes->GetAgglomerate(d)) continue; + if (reserved[c] || reserved[d]) continue; + + /*--- Geometrical quality check ---*/ + if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config) || + !GeometricalCheck(c, fine_grid, config) || !GeometricalCheck(d, fine_grid, config)) + continue; + + /*--- Guard against duplicate indices ---*/ + if (a == b || a == c || a == d || b == c || b == d || c == d) { + line_processed[li1] = line_processed[li2_best] = 1; + continue; + } + + /*--- Create 4-child coarse CV (isotropic agglomeration) ---*/ + fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 0, a); + fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 1, b); + fine_grid->nodes->SetParent_CV(c, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 2, c); + fine_grid->nodes->SetParent_CV(d, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 3, d); + nodes->SetnChildren_CV(Index_CoarseCV, 4); + + /*--- Debug output: show CV creation details ---*/ + if (DEBUG_OUTPUT && Index_CoarseCV < starting_Index_CoarseCV + DEBUG_CV_LIMIT) { + const auto* coord_a = fine_grid->nodes->GetCoord(a); + const auto* coord_b = fine_grid->nodes->GetCoord(b); + cout << " CV " << Index_CoarseCV << " (ISO): nodes " << a << "+" << b << "+" << c << "+" << d + << " | lines[" << li1 << "][" << idx1 << "," << idx2 << "]+lines[" << li2_best << "][" << idx1 << "," << idx2 << "]" + << " | coord_a=(" << coord_a[0] << "," << coord_a[1] << ")" + << " coord_b=(" << coord_b[0] << "," << coord_b[1] << ")" << endl; } + + reserved[a] = reserved[b] = reserved[c] = reserved[d] = 1; + MGQueue_InnerCV.RemoveCV(a); + MGQueue_InnerCV.RemoveCV(b); + MGQueue_InnerCV.RemoveCV(c); + MGQueue_InnerCV.RemoveCV(d); + + Index_CoarseCV++; + line_processed[li1] = line_processed[li2_best] = 1; + any_work = true; } } else { - /*--- ANISOTROPIC MODE: Pair nodes at SAME position on DIFFERENT lines ---*/ - for (auto& [parent, line_ids] : parent_to_lines) { - if (line_ids.size() < 2) continue; - - /*--- Pair consecutive lines at the same position ---*/ - for (size_t k = 0; k + 1 < line_ids.size(); k += 2) { - const auto li1 = line_ids[k]; - const auto li2 = line_ids[k + 1]; - if (line_processed[li1] || line_processed[li2]) continue; - - const auto& L1 = lines[li1]; - const auto& L2 = lines[li2]; - const auto pos = 1 + position_idx; - if (L1.size() <= pos || L2.size() <= pos) continue; - - const auto a = L1[pos]; - const auto b = L2[pos]; - - /*--- Skip if any node is already claimed ---*/ - if (fine_grid->nodes->GetAgglomerate(a) || fine_grid->nodes->GetAgglomerate(b)) continue; - if (reserved[a] || reserved[b]) continue; - - /*--- Geometrical quality check ---*/ - if (!GeometricalCheck(a, fine_grid, config) || !GeometricalCheck(b, fine_grid, config)) continue; - - /*--- Create 2-child coarse CV (anisotropic: same position, different lines) ---*/ - fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 0, a); - fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); - nodes->SetChildren_CV(Index_CoarseCV, 1, b); - nodes->SetnChildren_CV(Index_CoarseCV, 2); - - reserved[a] = reserved[b] = 1; - MGQueue_InnerCV.RemoveCV(a); - MGQueue_InnerCV.RemoveCV(b); - - Index_CoarseCV++; - line_processed[li1] = line_processed[li2] = 1; - any_work = true; + /*--- ANISOTROPIC MODE: Pair nodes at SAME position on DIFFERENT lines + Use spatial neighbor search to pair adjacent lines. ---*/ + for (auto li1 : active_lines) { + if (line_processed[li1]) continue; + + const auto& L1 = lines[li1]; + const auto pos = 1 + position_idx; + if (L1.size() <= pos) continue; + + const auto a = L1[pos]; + if (fine_grid->nodes->GetAgglomerate(a)) continue; + if (reserved[a]) continue; + if (!GeometricalCheck(a, fine_grid, config)) continue; + + /*--- Find nearest neighbor line by checking mesh neighbors of node 'a' ---*/ + unsigned long li2_best = std::numeric_limits::max(); + for (auto neighbor_point : fine_grid->nodes->GetPoints(a)) { + /*--- Check if this neighbor belongs to another unprocessed line at same position ---*/ + for (auto li2 : active_lines) { + if (li2 == li1 || line_processed[li2]) continue; + const auto& L2 = lines[li2]; + if (L2.size() <= pos) continue; + const auto b = L2[pos]; + if (b == neighbor_point) { + li2_best = li2; + break; + } + } + if (li2_best != std::numeric_limits::max()) break; + } + + if (li2_best == std::numeric_limits::max()) { + /*--- Debug: Line couldn't find a neighbor ---*/ + if (DEBUG_OUTPUT && position_idx < 3) { + cout << " Line " << li1 << " at pos=" << pos << " (node " << a << ") has NO neighbor line!" << endl; + } + continue; + } + + const auto& L2 = lines[li2_best]; + const auto b = L2[pos]; + + /*--- Skip if partner is already claimed ---*/ + if (fine_grid->nodes->GetAgglomerate(b)) continue; + if (reserved[b]) continue; + + /*--- Geometrical quality check ---*/ + if (!GeometricalCheck(b, fine_grid, config)) continue; + + /*--- Debug: Check line distance and neighbor relationships ---*/ + if (DEBUG_OUTPUT && Index_CoarseCV < starting_Index_CoarseCV + DEBUG_CV_LIMIT) { + /*--- Measure distance between wall vertices of the two lines ---*/ + const auto wall_a = lines[li1][0]; + const auto wall_b = lines[li2_best][0]; + const auto* coord_wall_a = fine_grid->nodes->GetCoord(wall_a); + const auto* coord_wall_b = fine_grid->nodes->GetCoord(wall_b); + su2double wall_dist = sqrt(pow(coord_wall_a[0] - coord_wall_b[0], 2) + + pow(coord_wall_a[1] - coord_wall_b[1], 2)); + + /*--- Check if wall vertices are neighbors ---*/ + bool walls_are_neighbors = false; + for (auto neighbor : fine_grid->nodes->GetPoints(wall_a)) { + if (neighbor == wall_b) { walls_are_neighbors = true; break; } + } + + cout << " Pairing lines " << li1 << " + " << li2_best << " at pos=" << pos + << " | wall_dist=" << wall_dist << " | walls_neighbors=" << (walls_are_neighbors ? "YES" : "NO") << endl; } + + /*--- Create 2-child coarse CV (anisotropic: same position, different lines) ---*/ + fine_grid->nodes->SetParent_CV(a, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 0, a); + fine_grid->nodes->SetParent_CV(b, Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, 1, b); + nodes->SetnChildren_CV(Index_CoarseCV, 2); + + /*--- Debug output: show CV creation details ---*/ + if (DEBUG_OUTPUT && Index_CoarseCV < starting_Index_CoarseCV + DEBUG_CV_LIMIT) { + const auto* coord_a = fine_grid->nodes->GetCoord(a); + const auto* coord_b = fine_grid->nodes->GetCoord(b); + su2double dist = sqrt(pow(coord_a[0] - coord_b[0], 2) + pow(coord_a[1] - coord_b[1], 2)); + bool are_neighbors = false; + for (auto neighbor : fine_grid->nodes->GetPoints(a)) { + if (neighbor == b) { are_neighbors = true; break; } + } + cout << " CV " << Index_CoarseCV << " (ANISO): nodes " << a << "+" << b + << " | lines[" << li1 << "][" << pos << "]+lines[" << li2_best << "][" << pos << "]" + << " | dist=" << dist << " | neighbors=" << (are_neighbors ? "YES" : "NO") + << " | coords A=(" << coord_a[0] << "," << coord_a[1] << ")" + << " B=(" << coord_b[0] << "," << coord_b[1] << ")" << endl; + } + + reserved[a] = reserved[b] = 1; + MGQueue_InnerCV.RemoveCV(a); + MGQueue_InnerCV.RemoveCV(b); + + Index_CoarseCV++; + line_processed[li1] = line_processed[li2_best] = 1; + any_work = true; } } @@ -1544,4 +1738,62 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, } if (!any_more) break; } + + /*--- Count how many CVs and nodes were created ---*/ + const auto nCVs_created = Index_CoarseCV - starting_Index_CoarseCV; + unsigned long nNodes_claimed = 0; + unsigned long nNodes_on_lines = 0; + unsigned long nNodes_unpaired = 0; + + for (const auto& L : lines) { + for (size_t i = 1; i < L.size(); ++i) { // Skip wall node at [0] + nNodes_on_lines++; + if (!reserved[L[i]]) nNodes_unpaired++; + } + } + + for (unsigned long i = 0; i < nPointFine; ++i) { + if (reserved[i]) nNodes_claimed++; + } + + if (rank == MASTER_NODE) { + cout << " Created " << nCVs_created << " coarse CVs from " << nNodes_claimed << " fine nodes." << endl; + cout << " Nodes on implicit lines: " << nNodes_on_lines << " (paired=" << (nNodes_on_lines - nNodes_unpaired) + << ", unpaired=" << nNodes_unpaired << ")" << endl; + + if (nNodes_unpaired > 0) { + cout << " WARNING: " << nNodes_unpaired << " nodes on implicit lines were left unpaired!" << endl; + cout << " These will be processed by domain agglomeration (may create wrong orientation)." << endl; + + /*--- Show first few unpaired nodes ---*/ + unsigned long count = 0; + for (size_t li = 0; li < lines.size() && count < 10; ++li) { + const auto& L = lines[li]; + for (size_t i = 1; i < L.size() && count < 10; ++i) { + if (!reserved[L[i]]) { + cout << " Unpaired: line " << li << " node " << L[i] << " at position " << i << endl; + count++; + } + } + } + } + if (ISOTROPIC) { + cout << " Expected ratio: ~4 nodes per CV (actual: " << std::setprecision(2) << std::fixed + << (nCVs_created > 0 ? su2double(nNodes_claimed) / su2double(nCVs_created) : 0.0) << ")" << endl; + } else { + cout << " Expected ratio: ~2 nodes per CV (actual: " << std::setprecision(2) << std::fixed + << (nCVs_created > 0 ? su2double(nNodes_claimed) / su2double(nCVs_created) : 0.0) << ")" << endl; + } + } + + /*--- Verify all claimed nodes are properly marked as agglomerated ---*/ + unsigned long mismatches = 0; + for (unsigned long i = 0; i < nPointFine; ++i) { + if (reserved[i] && !fine_grid->nodes->GetAgglomerate(i)) { + mismatches++; + } + } + if (mismatches > 0 && rank == MASTER_NODE) { + cout << " WARNING: " << mismatches << " nodes marked as reserved but not agglomerated!" << endl; + } } diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 1c0b35ad53d..cb19c8db88d 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -154,10 +154,11 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, /*--- Full MG: advance to the next finer grid after a fixed number of * outer iterations on the current coarsest active level. - * We use 100 iterations per level (nMGLevels levels total) ---*/ + * The number of iterations per level is controlled by MG_STARTUP_ITER config option. ---*/ + const unsigned long startup_iter = config[iZone]->GetMGOptions().MG_Startup_Iter; const bool Convergence_FullMG = FullMG && (FinestMesh != MESH_0) && - (config[iZone]->GetInnerIter() % 100 == 99); + (config[iZone]->GetInnerIter() % startup_iter == startup_iter - 1); if (!config[iZone]->GetRestart() && FullMG && direct && ( Convergence_FullMG && (FinestMesh != MESH_0 ))) { @@ -197,10 +198,10 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, passivedouble CFL_local = cfl_base; for (unsigned short iMesh = FinestMesh; iMesh < nMGLevels; ++iMesh) { const unsigned short lvl = iMesh + 1; - /*--- Use per-level scaling factor; clamp to (0,1] to prevent coarse CFL from - * exceeding the fine CFL. Index into cflScaling is iMesh (0-based transition). ---*/ + /*--- Use per-level scaling factor to increase coarse CFL (allows values > 1.0). + * Index into cflScaling is iMesh (0-based transition). ---*/ const passivedouble scale = (iMesh < cflScaling.size()) - ? max(passivedouble{1e-6}, min(passivedouble{1.0}, SU2_TYPE::GetValue(cflScaling[iMesh]))) + ? max(passivedouble{1e-6}, SU2_TYPE::GetValue(cflScaling[iMesh])) : passivedouble{0.25}; CFL_local *= scale; config[iZone]->SetCFL(lvl, CFL_local); diff --git a/config_template.cfg b/config_template.cfg index 9e2fe09703d..614338c5ba8 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1724,6 +1724,10 @@ MG_IMPLICIT_LINES_MAX_LENGTH= 20 % Anisotropic (NO): Pair cells normal to wall (2 cells per coarse CV, reduces mesh ~2x) % Isotropic (YES): Pair cells in all directions (4 cells per coarse CV, reduces mesh ~4x) MG_IMPLICIT_LINES_ISOTROPIC= NO +% +% Number of iterations on coarsest mesh during Full Multigrid (FMG) startup phase. +% After this many iterations, solution is prolongated to finer mesh (default 100). +MG_STARTUP_ITER= 100 % -------------------------- MESH SMOOTHING -----------------------------% % From 3b9e9dfbbb5e12d1196b0425c052cbf9b2f2bd8c Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 1 Aug 2026 16:45:31 +0200 Subject: [PATCH 3/7] fix multigrid turbulence --- Common/include/option_structure.hpp | 1 + Common/src/CConfig.cpp | 2 + .../integration/CMultiGridIntegration.hpp | 20 +++ .../src/integration/CMultiGridIntegration.cpp | 131 +++++++++++++++++- SU2_CFD/src/iteration/CFluidIteration.cpp | 15 +- SU2_CFD/src/solvers/CSolverFactory.cpp | 2 +- SU2_CFD/src/solvers/CTurbSASolver.cpp | 2 +- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 2 +- SU2_CFD/src/variables/CTurbVariable.cpp | 12 ++ config_template.cfg | 4 + 10 files changed, 183 insertions(+), 8 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index e277e514332..14e09e9158f 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1128,6 +1128,7 @@ struct CMGOptions { unsigned long MG_Implicit_Lines_MaxLength{20}; /*!< \brief Maximum nodes on a wall-normal implicit line (including wall seed). */ bool MG_Implicit_Lines_Isotropic{false}; /*!< \brief Use isotropic (vs anisotropic) agglomeration along implicit lines. */ unsigned long MG_Startup_Iter{100}; /*!< \brief Number of iterations on coarsest mesh during FMG startup phase. */ + bool TurbMG{false}; /*!< \brief Run turbulence equations through a FAS MG V-cycle instead of single-grid. */ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 2d68a229b27..208767d6db5 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2073,6 +2073,8 @@ void CConfig::SetConfig_Options() { addBoolOption("MG_IMPLICIT_LINES_ISOTROPIC", MGOptions.MG_Implicit_Lines_Isotropic, false); /*!\brief MG_STARTUP_ITER\n DESCRIPTION: Number of iterations on the coarsest mesh during Full Multigrid (FMG) startup phase before advancing to finer meshes. DEFAULT: 100 \ingroup Config*/ addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100); + /*!\brief MG_TURB\n DESCRIPTION: Run turbulence equations through a FAS Multigrid V-cycle instead of single-grid. DEFAULT: NO \ingroup Config*/ + addBoolOption("MG_TURB", MGOptions.TurbMG, false); /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i). If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/ addDoubleListOption("MG_CFL_SCALING", nMG_CflScaling_p, MG_CflScaling_p); diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp index 0672d1ae5ff..37f68c33f8f 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -248,6 +248,26 @@ class CMultiGridIntegration final : public CIntegration { passivedouble lastRMS[2], char& exitReason, passivedouble& worstStepRatio, unsigned short& worstStep); + /*! + * \brief Restrict turbulent eddy viscosity from fine to coarser mesh levels. + * + * After a turbulence FAS V-cycle completes, this function volume-weights restricts + * mu_t from the finest mesh down to all coarser levels. The flow solver on the next + * outer iteration uses these restricted mu_t values at every coarse level for the + * eddy-viscosity coupling. This ensures consistency between flow and turbulence + * solutions across the multigrid hierarchy. + * + * \param[in] geometry - Geometry hierarchy for one zone/instance (all levels). + * \param[in] solver - Solver hierarchy for one zone/instance (all levels). + * \param[in] config - Problem configuration. + * \param[in] FinestMesh - Current finest active mesh index. + * \param[in] nMGLevels - Total number of MG levels. + */ + void RestrictTurbEddyViscToCoarseLevels(CGeometry** geometry, CSolver*** solver, + CConfig* config, + unsigned short FinestMesh, + unsigned short nMGLevels); + static constexpr int MAX_MG_LEVELS = 10; /*--- Early-exit smoothing state (shared across OMP threads via master write + barrier). ---*/ diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index cb19c8db88d..6ebbfc9eeb5 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -28,6 +28,11 @@ #include "../../include/integration/CMultiGridIntegration.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" +#include +#include +#include + +using namespace std; namespace { @@ -160,7 +165,8 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, FullMG && (FinestMesh != MESH_0) && (config[iZone]->GetInnerIter() % startup_iter == startup_iter - 1); - if (!config[iZone]->GetRestart() && FullMG && direct && ( Convergence_FullMG && (FinestMesh != MESH_0 ))) { + if (!config[iZone]->GetRestart() && FullMG && direct && ( Convergence_FullMG && (FinestMesh != MESH_0 )) && + RunTime_EqSystem == RUNTIME_FLOW_SYS) { SetProlongated_Solution(RunTime_EqSystem, solver_container[iZone][iInst][FinestMesh-1][Solver_Position], @@ -169,6 +175,46 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, geometry[iZone][iInst][FinestMesh], config[iZone]); + /*--- Prolongate scalar solutions to the new finest mesh. + * All scalar solvers (turb, species, transition) run via SingleGrid_Iteration on + * GetFinestMesh(). Only turbulence additionally restricts its field downward to + * coarser meshes; no scalar ever propagates upward to finer meshes. Consequently + * meshes finer than FinestMesh hold their iter-0 startup values for the entire + * warmup phase. When FinestMesh is decremented these stale fields cause a large + * transient (e.g. +3 decade regression in rms[nu]). Prolongating here mirrors + * what SetProlongated_Solution does for the flow and eliminates the regression. ---*/ + if (config[iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE) { + SetProlongated_Solution(RUNTIME_TURB_SYS, + solver_container[iZone][iInst][FinestMesh-1][TURB_SOL], + solver_container[iZone][iInst][FinestMesh][TURB_SOL], + geometry[iZone][iInst][FinestMesh-1], + geometry[iZone][iInst][FinestMesh], + config[iZone]); + /*--- Recompute mu_t on the new finest mesh from the prolongated nu_tilde/k/omega. ---*/ + solver_container[iZone][iInst][FinestMesh-1][TURB_SOL]->Postprocessing( + geometry[iZone][iInst][FinestMesh-1], + solver_container[iZone][iInst][FinestMesh-1], + config[iZone], FinestMesh-1); + } + + if (config[iZone]->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { + SetProlongated_Solution(RUNTIME_TRANS_SYS, + solver_container[iZone][iInst][FinestMesh-1][TRANS_SOL], + solver_container[iZone][iInst][FinestMesh][TRANS_SOL], + geometry[iZone][iInst][FinestMesh-1], + geometry[iZone][iInst][FinestMesh], + config[iZone]); + } + + if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) { + SetProlongated_Solution(RUNTIME_SPECIES_SYS, + solver_container[iZone][iInst][FinestMesh-1][SPECIES_SOL], + solver_container[iZone][iInst][FinestMesh][SPECIES_SOL], + geometry[iZone][iInst][FinestMesh-1], + geometry[iZone][iInst][FinestMesh], + config[iZone]); + } + SU2_OMP_SAFE_GLOBAL_ACCESS(config[iZone]->SubtractFinestMesh();) } @@ -176,11 +222,45 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, FinestMesh = config[iZone]->GetFinestMesh(); + /*--- For turbulence MG: before descending to coarse levels, ensure mu_t is computed + * at the finest level and restricted to all coarser levels. This prevents inf + * residuals from coarse-level turbulence solves using stale/uninitialized mu_t. ---*/ + if (RunTime_EqSystem == RUNTIME_TURB_SYS && + config[iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE) { + + solver_container[iZone][iInst][FinestMesh][TURB_SOL]->Postprocessing( + geometry[iZone][iInst][FinestMesh], + solver_container[iZone][iInst][FinestMesh], + config[iZone], FinestMesh); + + RestrictTurbEddyViscToCoarseLevels(geometry[iZone][iInst], + solver_container[iZone][iInst], + config[iZone], FinestMesh, + config[iZone]->GetnMGLevels()); + } + /*--- Perform the Full Approximation Scheme multigrid ---*/ MultiGrid_Cycle(geometry, solver_container, numerics_container, config, FinestMesh, RecursiveParam, RunTime_EqSystem, iZone, iInst); + /*--- After a turb FAS V-cycle: recompute mu_t at the finest active level from the updated + * nu_hat/k/omega and restrict it to all coarser levels. The flow FAS on the NEXT outer + * iteration uses these mu_t values at every coarse level for the eddy-viscosity coupling. + * (Postprocessing was already called on FinestMesh inside the last PreSmoothing step of + * MultiGrid_Cycle; we call it once more to be safe after the V-cycle correction is applied.) ---*/ + if (RunTime_EqSystem == RUNTIME_TURB_SYS && + config[iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE) { + solver_container[iZone][iInst][FinestMesh][TURB_SOL]->Postprocessing( + geometry[iZone][iInst][FinestMesh], + solver_container[iZone][iInst][FinestMesh], + config[iZone], FinestMesh); + RestrictTurbEddyViscToCoarseLevels(geometry[iZone][iInst], + solver_container[iZone][iInst], + config[iZone], FinestMesh, + config[iZone]->GetnMGLevels()); + } + /*--- Adapt coarse-grid CFL once per cycle using smoothing residuals gathered during the cycle. ---*/ const unsigned short nMGLevels = config[iZone]->GetnMGLevels(); BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS @@ -284,8 +364,13 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, return ss.str(); }; + const string eqName = (RunTime_EqSystem == RUNTIME_FLOW_SYS) ? "Flow" : + (RunTime_EqSystem == RUNTIME_TURB_SYS) ? "Turb" : + (RunTime_EqSystem == RUNTIME_SPECIES_SYS) ? "Species" : + (RunTime_EqSystem == RUNTIME_TRANS_SYS) ? "Trans" : "Other"; + PrintingToolbox::CTablePrinter table(&std::cout); - table.AddColumn("Smoother", 13); + table.AddColumn("Smoother [" + eqName + "]", 13 + 7); for (unsigned short i = 0; i <= nMGLevels; ++i) table.AddColumn("Level " + std::to_string(i), 38); table.PrintHeader(); @@ -412,6 +497,15 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, solver_coarse->Preprocessing(geometry_coarse, solver_container_coarse, config, iMesh+1, NO_RK_ITER, RunTime_EqSystem, false); + /*--- For turbulence: ensure flow primitives (density, laminar viscosity) are updated on the + * coarse level from the restricted conservative variables, THEN compute mu_t from the + * newly restricted turbulence variables. This ensures turbulence Postprocessing reads + * valid flow data and Space_Integration uses correct eddy viscosity. ---*/ + if (RunTime_EqSystem == RUNTIME_TURB_SYS && config->GetKind_Turb_Model() != TURB_MODEL::NONE) { + solver_container_coarse[FLOW_SOL]->Preprocessing(geometry_coarse, solver_container_coarse, config, iMesh+1, NO_RK_ITER, RUNTIME_FLOW_SYS, false); + solver_coarse->Postprocessing(geometry_coarse, solver_container_coarse, config, iMesh+1); + } + Space_Integration(geometry_coarse, solver_container_coarse, numerics_coarse, config, iMesh+1, NO_RK_ITER, RunTime_EqSystem); /*--- Compute $P_(k+1) = I^(k+1)_k(r_k) - r_(k+1) ---*/ @@ -1095,3 +1189,36 @@ void CMultiGridIntegration::Adjoint_Setup(CGeometry ****geometry, CSolver *****s } } + +void CMultiGridIntegration::RestrictTurbEddyViscToCoarseLevels(CGeometry** geometry, CSolver*** solver, + CConfig* config, + unsigned short FinestMesh, + unsigned short nMGLevels) { + SU2_ZONE_SCOPED + + for (unsigned short iMesh = FinestMesh; iMesh < nMGLevels; iMesh++) { + + CGeometry* geo_fine = geometry[iMesh]; + CGeometry* geo_coarse = geometry[iMesh + 1]; + CSolver* sol_fine = solver[iMesh][TURB_SOL]; + CSolver* sol_coarse = solver[iMesh + 1][TURB_SOL]; + + /*--- Volume-weighted restriction of mu_t from fine to coarse. ---*/ + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + + const su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); + su2double EddyVisc = 0.0; + + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + su2double Area_Children = geo_fine->nodes->GetVolume(Point_Fine); + su2double mu_t = sol_fine->GetNodes()->GetmuT(Point_Fine); + EddyVisc += mu_t * Area_Children / Area_Parent; + } + + sol_coarse->GetNodes()->SetmuT(Point_Coarse, EddyVisc); + } + END_SU2_OMP_FOR + } +} diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index a79a9f1f8f4..bdfd71fb8fc 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -82,7 +82,7 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe integration[val_iZone][val_iInst][FLOW_SOL]->MultiGrid_Iteration(geometry, solver, numerics, config, RUNTIME_FLOW_SYS, val_iZone, val_iInst); - /*--- If the flow integration is not fully coupled, run the various single grid integrations. ---*/ + /*--- If the flow integration is not fully coupled, run the various single/multi-grid integrations. ---*/ if (config[val_iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE && !frozen_visc) { @@ -95,14 +95,22 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe } /*--- Solve the turbulence model ---*/ + /*--- Use multigrid if MG_TURB is enabled, otherwise use single-grid. ---*/ config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_TURB_SYS); - integration[val_iZone][val_iInst][TURB_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, - RUNTIME_TURB_SYS, val_iZone, val_iInst); + + if (config[val_iZone]->GetMGOptions().TurbMG) { + integration[val_iZone][val_iInst][TURB_SOL]->MultiGrid_Iteration(geometry, solver, numerics, config, + RUNTIME_TURB_SYS, val_iZone, val_iInst); + } else { + integration[val_iZone][val_iInst][TURB_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, + RUNTIME_TURB_SYS, val_iZone, val_iInst); + } } if (config[val_iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) { config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_SPECIES_SYS); + /*--- Use multigrid if MG_SPECIES is enabled (future feature), otherwise use single-grid. ---*/ integration[val_iZone][val_iInst][SPECIES_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_SPECIES_SYS, val_iZone, val_iInst); @@ -119,6 +127,7 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe if (config[val_iZone]->GetWeakly_Coupled_Heat()) { config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_HEAT_SYS); + /*--- Use multigrid if MG_HEAT is enabled (future feature), otherwise use single-grid. ---*/ integration[val_iZone][val_iInst][HEAT_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_HEAT_SYS, val_iZone, val_iInst); } diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index ee798c1384b..00e7064d31c 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -313,7 +313,7 @@ CSolver* CSolverFactory::CreateSubSolver(SUB_SOLVER_TYPE kindSolver, CSolver **s case SUB_SOLVER_TYPE::TURB_SA: case SUB_SOLVER_TYPE::TURB_SST: genericSolver = CreateTurbSolver(kindTurbModel, solver, geometry, config, iMGLevel, false); - metaData.integrationType = INTEGRATION_TYPE::SINGLEGRID; + metaData.integrationType = config->GetMGOptions().TurbMG ? INTEGRATION_TYPE::MULTIGRID : INTEGRATION_TYPE::SINGLEGRID; break; case SUB_SOLVER_TYPE::TEMPLATE: genericSolver = new CTemplateSolver(geometry, config); diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index b58afe4c9b3..714d1972695 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -65,7 +65,7 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor /*--- Single grid simulation ---*/ - if (iMesh == MESH_0 || config->GetMGCycle() == MG_CYCLE::FULL) { + if (iMesh == MESH_0 || config->GetMGCycle() == MG_CYCLE::FULL || config->GetMGOptions().TurbMG) { /*--- Define some auxiliar vector related with the residual ---*/ diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 579b5f3cf72..97b4ee1d78c 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -59,7 +59,7 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- Single grid simulation ---*/ - if (iMesh == MESH_0 || config->GetMGCycle() == MG_CYCLE::FULL) { + if (iMesh == MESH_0 || config->GetMGCycle() == MG_CYCLE::FULL || config->GetMGOptions().TurbMG) { /*--- Define some auxiliary vector related with the residual ---*/ diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp index 139223055a4..62b87f9363d 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -35,6 +35,18 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned turb_index.resize(nPoint) = su2double(1.0); intermittency.resize(nPoint) = su2double(1.0); + /*--- Allocate residual structures for multigrid (required for turbulence MG). ---*/ + Res_TruncError.resize(nPoint, nVar) = su2double(0.0); + + /*--- Allocate smoothing arrays if correction smoothing is enabled at any MG level. ---*/ + for (unsigned long iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) { + if (config->GetMGOptions().MG_CorrecSmooth[iMesh] > 0) { + Residual_Sum.resize(nPoint, nVar); + Residual_Old.resize(nPoint, nVar); + break; + } + } + } void CTurbVariable::RegisterEddyViscosity(bool input) { diff --git a/config_template.cfg b/config_template.cfg index 45a822755c9..2c7d882c11f 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1732,6 +1732,10 @@ MG_IMPLICIT_LINES_ISOTROPIC= NO % Number of iterations on coarsest mesh during Full Multigrid (FMG) startup phase. % After this many iterations, solution is prolongated to finer mesh (default 100). MG_STARTUP_ITER= 100 +% +% Run turbulence equations through a FAS Multigrid V-cycle (YES, NO) +% When disabled, turbulence is solved with single-grid only (default NO). +MG_TURB= NO % -------------------------- MESH SMOOTHING -----------------------------% % From 7e1f670ab24b22f3b1e3a9d2cccce5d4ad531d39 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 2 Aug 2026 13:47:34 +0200 Subject: [PATCH 4/7] small flow solver update --- .../integration/CMultiGridIntegration.hpp | 26 +++ .../src/integration/CMultiGridIntegration.cpp | 180 ++++++++++++++++-- 2 files changed, 190 insertions(+), 16 deletions(-) diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp index 37f68c33f8f..1b5c7476bf6 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -81,6 +81,18 @@ class CMultiGridIntegration final : public CIntegration { void SetForcing_Term(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, unsigned short iMesh); + /*! + * \brief Restrict the fine-grid residual defect to the coarse-grid FAS forcing term. + * \param[in] sol_fine - Pointer to the solution on the fine grid. + * \param[in] sol_coarse - Pointer to the solution on the coarse grid. + * \param[in] geo_fine - Geometrical definition of the fine grid. + * \param[in] geo_coarse - Geometrical definition of the coarse grid. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void RestrictResidualToCoarseGrid(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, + CGeometry *geo_coarse, CConfig *config, unsigned short iMesh); + /*! * \brief Add the truncation error to the residual. * \param[in] geometry - Geometrical definition of the problem. @@ -180,6 +192,20 @@ class CMultiGridIntegration final : public CIntegration { void GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + /*! + * \brief Prolongate the coarse-grid state correction back to the fine-grid residual correction. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] sol_fine - Pointer to the solution on the fine grid. + * \param[in] sol_coarse - Pointer to the solution on the coarse grid. + * \param[in] geo_fine - Geometrical definition of the fine grid. + * \param[in] geo_coarse - Geometrical definition of the coarse grid. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void ProlongateCorrectionToFineGrid(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, + unsigned short iMesh); + /*! * \brief Do an implicit smoothing of the prolongated correction. * \param[in] RunTime_EqSystem - System of equations which is going to be solved. diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 6ebbfc9eeb5..91097a635d2 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -460,7 +460,8 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, PreSmoothing(RunTime_EqSystem, geometry, solver_container, config_container, solver_fine, numerics_fine, geometry_fine, solver_container_fine, config, iMesh, iZone, iRKLimit); - /*--- Compute Forcing Term $P_(k+1) = I^(k+1)_k(P_k+F_k(u_k))-F_(k+1)(I^(k+1)_k u_k)$ and update solution for multigrid ---*/ + /*--- Assemble the coarse-grid FAS defect term by restricting the fine-grid residual defect, + * solving the coarse-grid state, and prolongating only the state correction back to the fine grid. ---*/ if ( iMesh < config->GetnMGLevels() ) { @@ -483,15 +484,16 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, Space_Integration(geometry_fine, solver_container_fine, numerics_fine, config, iMesh, NO_RK_ITER, RunTime_EqSystem); - /*--- LinSysRes = R(u_N) here, before tau is added by SetResidual_Term. ---*/ + /*--- LinSysRes = R(u_N) here, before the fine-grid defect term is assembled. ---*/ BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { lastPreSmoothRMS[iMesh][1] = ComputeLinSysResRMS(solver_fine); } END_SU2_OMP_SAFE_GLOBAL_ACCESS + /*--- Assemble the fine-grid defect term that will be restricted to the coarse-grid FAS problem. ---*/ SetResidual_Term(geometry_fine, solver_fine); - /*--- Compute $r_(k+1) = F_(k+1)(I^(k+1)_k u_k)$ ---*/ + /*--- Restrict the fine-grid state to the coarse grid and initialize the coarse-grid state. ---*/ SetRestricted_Solution(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); @@ -508,9 +510,12 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, Space_Integration(geometry_coarse, solver_container_coarse, numerics_coarse, config, iMesh+1, NO_RK_ITER, RunTime_EqSystem); - /*--- Compute $P_(k+1) = I^(k+1)_k(r_k) - r_(k+1) ---*/ - - SetForcing_Term(solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, iMesh+1); + /*--- Restrict the fine-grid residual defect to the coarse-grid FAS forcing term. ---*/ + if (RunTime_EqSystem == RUNTIME_FLOW_SYS) { + RestrictResidualToCoarseGrid(solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, iMesh+1); + } else { + SetForcing_Term(solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, iMesh+1); + } /*--- Restore the time integration settings. ---*/ @@ -532,9 +537,12 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, iMesh+1, nextRecurseParam, RunTime_EqSystem, iZone, iInst); } - /*--- Compute prolongated solution, and smooth the correction $u^(new)_k = u_k + Smooth(I^k_(k+1)(u_(k+1)-I^(k+1)_k u_k))$ ---*/ - - GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); + /*--- Compute the coarse-grid state correction and prolongate it back to the fine grid. ---*/ + if (RunTime_EqSystem == RUNTIME_FLOW_SYS) { + ProlongateCorrectionToFineGrid(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, iMesh); + } else { + GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); + } const auto& mgOpts = config->GetMGOptions(); SmoothProlongated_Correction(RunTime_EqSystem, solver_fine, geometry_fine, mgOpts.MG_CorrecSmooth[iMesh], mgOpts.MG_Smooth_Coeff, config, iMesh); @@ -823,6 +831,8 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ if (val_nSmooth == 0) return; const unsigned short nVar = solver->GetnVar(); + const bool use_conservative_damping = (nVar <= 2); + const su2double turbulence_base_damping = 0.50; SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads())) for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); iPoint++) { @@ -862,8 +872,11 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ const auto* Residual_Sum = solver->GetNodes()->GetResidual_Sum(iPoint); const auto* Residual_Old = solver->GetNodes()->GetResidual_Old(iPoint); - for (auto iVar = 0u; iVar < nVar; iVar++) - solver->LinSysRes(iPoint,iVar) = (Residual_Old[iVar] + val_smooth_coeff*Residual_Sum[iVar])*factor; + for (auto iVar = 0u; iVar < nVar; iVar++) { + su2double smoothed = (Residual_Old[iVar] + val_smooth_coeff*Residual_Sum[iVar])*factor; + if (use_conservative_damping) smoothed *= turbulence_base_damping; + solver->LinSysRes(iPoint,iVar) = smoothed; + } } END_SU2_OMP_FOR @@ -890,6 +903,11 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ if (config->GetMGOptions().MG_Smooth_Output) { const su2double res = sqrt(solver->LinSysRes.squaredNorm() / (nVar * geometry->GetGlobal_nPointDomain())); SU2_OMP_SAFE_GLOBAL_ACCESS(lastCorrecSmoothRMS[iMesh][1] = SU2_TYPE::GetValue(res);) + + if (SU2_MPI::GetRank() == MASTER_NODE && use_conservative_damping) { + cout << "[MG CORR-SMOOTH] turbulence nSmooth=" << val_nSmooth + << " norm=" << res << "\n"; + } } } @@ -898,25 +916,129 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet SU2_ZONE_SCOPED const unsigned short nVar = sol_fine->GetnVar(); + const bool use_conservative_damping = (nVar <= 2); + const su2double base_damping = use_conservative_damping ? 0.50 : 1.0; + const su2double wall_damping = use_conservative_damping ? 0.25 : 1.0; /*--- Use the adaptive damping factor uniformly across all prolongation levels. ---*/ const su2double factor = config->GetDamp_Correc_Prolong(); + vector isWall(geo_fine->GetnPoint(), false); + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) + if (config->GetViscous_Wall(iMarker)) + for (auto iVertex = 0ul; iVertex < geo_fine->nVertex[iMarker]; iVertex++) + isWall[geo_fine->vertex[iMarker][iVertex]->GetNode()] = true; + SU2_OMP_FOR_STAT(roundUpDiv(geo_fine->GetnPointDomain(), omp_get_num_threads())) for (auto Point_Fine = 0ul; Point_Fine < geo_fine->GetnPointDomain(); Point_Fine++) { auto* Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine); auto* Solution_Fine = sol_fine->GetNodes()->GetSolution(Point_Fine); + + su2double residualMag = 0.0; + su2double correctionMag = 0.0; for (auto iVar = 0u; iVar < nVar; iVar++) { /*--- Prevent a fine grid divergence due to a coarse grid divergence ---*/ - if (Residual_Fine[iVar] != Residual_Fine[iVar]) + if (Residual_Fine[iVar] != Residual_Fine[iVar]) { Residual_Fine[iVar] = 0.0; + } + const su2double corr = factor * Residual_Fine[iVar]; + residualMag = max(residualMag, fabs(Residual_Fine[iVar])); + correctionMag = max(correctionMag, fabs(corr)); + } + + su2double correctionScale = 1.0; + if (residualMag > 1e-30 && correctionMag > 1e-30) { + const su2double ratio = correctionMag / residualMag; + if (ratio > 2.0) { + correctionScale = 2.0 / ratio; + } + } + + const su2double localDamping = use_conservative_damping ? (isWall[Point_Fine] ? wall_damping : base_damping) : 1.0; + for (auto iVar = 0u; iVar < nVar; iVar++) { su2double correction = factor * Residual_Fine[iVar]; + correction *= localDamping; + correction *= correctionScale; + + if (!std::isfinite(correction)) { + correction = 0.0; + } + Solution_Fine[iVar] += correction; } } END_SU2_OMP_FOR + /*--- DIAGNOSTIC: log the max applied correction (factor * LinSysRes) at fine-grid wall points + * vs interior. ---*/ + if (config->GetMGOptions().MG_Smooth_Output && SU2_MPI::GetRank() == MASTER_NODE) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { + if (nVar > 2) { + su2double maxWall0 = 0.0, maxWallN = 0.0, maxWallMom = 0.0; + su2double maxInter0 = 0.0, maxInterN = 0.0, maxInterMom = 0.0; + su2double maxWallApply0 = 0.0, maxWallApplyN = 0.0, maxWallApplyMom = 0.0; + su2double maxInterApply0 = 0.0, maxInterApplyN = 0.0, maxInterApplyMom = 0.0; + + for (auto iPoint = 0ul; iPoint < geo_fine->GetnPointDomain(); iPoint++) { + const auto* corr = sol_fine->LinSysRes.GetBlock(iPoint); + const su2double localDamping = isWall[iPoint] ? wall_damping : base_damping; + const su2double applied0 = fabs(localDamping * factor * corr[0]); + const su2double appliedN = fabs(localDamping * factor * corr[nVar-1]); + su2double appliedMom = 0.0; + for (auto iVar = 1u; iVar < static_cast(nVar-1); iVar++) { + appliedMom = max(appliedMom, fabs(localDamping * factor * corr[iVar])); + } + + if (isWall[iPoint]) { + maxWall0 = max(maxWall0, fabs(factor * corr[0])); + maxWallN = max(maxWallN, fabs(factor * corr[nVar-1])); + maxWallMom = max(maxWallMom, fabs(factor * corr[0])); + maxWallApply0 = max(maxWallApply0, applied0); + maxWallApplyN = max(maxWallApplyN, appliedN); + maxWallApplyMom = max(maxWallApplyMom, appliedMom); + } else { + maxInter0 = max(maxInter0, fabs(factor * corr[0])); + maxInterN = max(maxInterN, fabs(factor * corr[nVar-1])); + maxInterMom = max(maxInterMom, fabs(factor * corr[0])); + maxInterApply0 = max(maxInterApply0, applied0); + maxInterApplyN = max(maxInterApplyN, appliedN); + maxInterApplyMom = max(maxInterApplyMom, appliedMom); + } + } + auto ratio = [](su2double w, su2double i) { return (i > 1e-30) ? w/i : 0.0; }; + cout << "[MG APPL wall/inter] rho=" << ratio(maxWallApply0, maxInterApply0) + << " mom=" << ratio(maxWallApplyMom, maxInterApplyMom) + << " E=" << ratio(maxWallApplyN, maxInterApplyN) + << " (raw wall/inter: rho=" << maxWall0 << "/" << maxInter0 + << ", E=" << maxWallN << "/" << maxInterN + << "; applied wall/inter: rho=" << maxWallApply0 << "/" << maxInterApply0 + << ", E=" << maxWallApplyN << "/" << maxInterApplyN << ")\n"; + } else { + su2double maxWall = 0.0, maxInter = 0.0; + su2double maxWallApply = 0.0, maxInterApply = 0.0; + for (auto iPoint = 0ul; iPoint < geo_fine->GetnPointDomain(); iPoint++) { + const auto* corr = sol_fine->LinSysRes.GetBlock(iPoint); + const su2double localDamping = isWall[iPoint] ? wall_damping : base_damping; + const su2double mag = fabs(factor * corr[0]); + const su2double appliedMag = fabs(localDamping * factor * corr[0]); + if (isWall[iPoint]) { + maxWall = max(maxWall, mag); + maxWallApply = max(maxWallApply, appliedMag); + } else { + maxInter = max(maxInter, mag); + maxInterApply = max(maxInterApply, appliedMag); + } + } + cout << "[MG TURB APPLY] damp(wall/interior)= " << wall_damping << "/" << base_damping + << " raw max(wall/interior)= " << maxWall << "/" << maxInter + << " applied max(wall/interior)= " << maxWallApply << "/" << maxInterApply << "\n"; + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + /*--- MPI the new interpolated solution ---*/ sol_fine->InitiateComms(geo_fine, config, MPI_QUANTITIES::SOLUTION); @@ -947,21 +1069,21 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar const unsigned short nVar = sol_coarse->GetnVar(); const su2double factor = config->GetDamp_Res_Restric(); - su2activevector Residual(nVar); + su2activevector RestrictedDefect(nVar); SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { sol_coarse->GetNodes()->SetRes_TruncErrorZero(Point_Coarse); - Residual = su2double(0); + RestrictedDefect = su2double(0); for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine); for (auto iVar = 0u; iVar < nVar; iVar++) - Residual[iVar] += factor * Residual_Fine[iVar]; + RestrictedDefect[iVar] += factor * Residual_Fine[iVar]; } - sol_coarse->GetNodes()->AddRes_TruncError(Point_Coarse, Residual.data()); + sol_coarse->GetNodes()->AddRes_TruncError(Point_Coarse, RestrictedDefect.data()); } END_SU2_OMP_FOR @@ -996,6 +1118,32 @@ void CMultiGridIntegration::SetResidual_Term(CGeometry *geometry, CSolver *solve } +void CMultiGridIntegration::RestrictResidualToCoarseGrid(CSolver *sol_fine, CSolver *sol_coarse, + CGeometry *geo_fine, CGeometry *geo_coarse, + CConfig *config, unsigned short iMesh) { + SU2_ZONE_SCOPED + + /*--- This is the standard FAS restriction step: the fine-grid defect is passed to the + * coarse-grid problem as a forcing term. The existing SetForcing_Term routine already + * implements the conservative volume-weighted transfer and the damping factor in the + * same way the original MG cycle expects. ---*/ + SetForcing_Term(sol_fine, sol_coarse, geo_fine, geo_coarse, config, iMesh); +} + +void CMultiGridIntegration::ProlongateCorrectionToFineGrid(unsigned short RunTime_EqSystem, CSolver *sol_fine, + CSolver *sol_coarse, CGeometry *geo_fine, + CGeometry *geo_coarse, CConfig *config, + unsigned short iMesh) { + SU2_ZONE_SCOPED + + /*--- This is the standard FAS prolongation step: build the coarse-grid state correction, + * then transfer that correction to the fine-grid residual correction. The original + * GetProlongated_Correction routine already performs this transfer in the correct form; + * the additional scaling here would be equivalent to changing the correction operator. + * Keep the transfer operator unchanged and let the existing damping path control the size. ---*/ + GetProlongated_Correction(RunTime_EqSystem, sol_fine, sol_coarse, geo_fine, geo_coarse, config); +} + void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { SU2_ZONE_SCOPED From 33f2eaf14d55ffea6e3ace0204ba006b73629a95 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 2 Aug 2026 14:45:17 +0200 Subject: [PATCH 5/7] small flow solver update --- .../src/integration/CMultiGridIntegration.cpp | 29 ++++++++++++++----- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 91097a635d2..12997293c87 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -58,6 +58,15 @@ static su2double applyGlobalTrend(su2double factor, passivedouble crossCycleRati return max(su2double{CLAMP_MIN}, min(su2double{CLAMP_MAX}, factor)); } +static su2double GetMGLevelCorrectionScale(unsigned short iMesh) { + switch (iMesh) { + case 0: return 1.00; + case 1: return 0.75; + case 2: return 0.50; + default: return 0.35; + } +} + inline passivedouble ComputeLinSysResRMS(const CSolver* solver) { passivedouble result = 0; for (unsigned short iVar = 0; iVar < solver->GetnVar(); ++iVar) { @@ -917,8 +926,9 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet const unsigned short nVar = sol_fine->GetnVar(); const bool use_conservative_damping = (nVar <= 2); - const su2double base_damping = use_conservative_damping ? 0.50 : 1.0; - const su2double wall_damping = use_conservative_damping ? 0.25 : 1.0; + const su2double levelScale = GetMGLevelCorrectionScale(iMesh); + const su2double base_damping = use_conservative_damping ? max(su2double{0.15}, 0.50 * levelScale) : 1.0; + const su2double wall_damping = use_conservative_damping ? max(su2double{0.10}, 0.25 * levelScale) : 1.0; /*--- Use the adaptive damping factor uniformly across all prolongation levels. ---*/ const su2double factor = config->GetDamp_Correc_Prolong(); @@ -948,10 +958,11 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet } su2double correctionScale = 1.0; + constexpr su2double maxAllowedRatio = 1.25; if (residualMag > 1e-30 && correctionMag > 1e-30) { const su2double ratio = correctionMag / residualMag; - if (ratio > 2.0) { - correctionScale = 2.0 / ratio; + if (ratio > maxAllowedRatio) { + correctionScale = maxAllowedRatio / ratio; } } @@ -1008,13 +1019,14 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet } } auto ratio = [](su2double w, su2double i) { return (i > 1e-30) ? w/i : 0.0; }; - cout << "[MG APPL wall/inter] rho=" << ratio(maxWallApply0, maxInterApply0) + cout << "[MG APPL L" << iMesh << " wall/inter] rho=" << ratio(maxWallApply0, maxInterApply0) << " mom=" << ratio(maxWallApplyMom, maxInterApplyMom) << " E=" << ratio(maxWallApplyN, maxInterApplyN) << " (raw wall/inter: rho=" << maxWall0 << "/" << maxInter0 << ", E=" << maxWallN << "/" << maxInterN << "; applied wall/inter: rho=" << maxWallApply0 << "/" << maxInterApply0 - << ", E=" << maxWallApplyN << "/" << maxInterApplyN << ")\n"; + << ", E=" << maxWallApplyN << "/" << maxInterApplyN + << "; levelScale=" << levelScale << ", damp(wall/inter)=" << wall_damping << "/" << base_damping << ")\n"; } else { su2double maxWall = 0.0, maxInter = 0.0; su2double maxWallApply = 0.0, maxInterApply = 0.0; @@ -1031,9 +1043,10 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet maxInterApply = max(maxInterApply, appliedMag); } } - cout << "[MG TURB APPLY] damp(wall/interior)= " << wall_damping << "/" << base_damping + cout << "[MG TURB APPLY L" << iMesh << "] damp(wall/interior)= " << wall_damping << "/" << base_damping << " raw max(wall/interior)= " << maxWall << "/" << maxInter - << " applied max(wall/interior)= " << maxWallApply << "/" << maxInterApply << "\n"; + << " applied max(wall/interior)= " << maxWallApply << "/" << maxInterApply + << " levelScale=" << levelScale << "\n"; } } END_SU2_OMP_SAFE_GLOBAL_ACCESS From f66d614becf4db5451bb6a6b5c5452807d83a83d Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 2 Aug 2026 19:06:22 +0200 Subject: [PATCH 6/7] cfl adaptation for full multigrid --- SU2_CFD/include/solvers/CSolver.hpp | 10 + .../src/integration/CMultiGridIntegration.cpp | 178 ++++++++++++++---- SU2_CFD/src/iteration/CFluidIteration.cpp | 9 +- 3 files changed, 151 insertions(+), 46 deletions(-) diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 942e2e25877..0b13f958867 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -389,6 +389,16 @@ class CSolver { */ inline su2double GetAvg_CFL_Local(void) const { return Avg_CFL_Local; } + /*!\ + * \brief Set min/max/avg local CFL summary statistics. + * \param[in] val_cfl - Uniform CFL value to report. + */ + inline void SetCFL_Local_Stats(su2double val_cfl) { + Min_CFL_Local = val_cfl; + Max_CFL_Local = val_cfl; + Avg_CFL_Local = val_cfl; + } + /*! * \brief Get the number of variables of the problem. */ diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 12997293c87..2429a0e33d1 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -31,6 +31,7 @@ #include #include #include +#include using namespace std; @@ -75,6 +76,96 @@ inline passivedouble ComputeLinSysResRMS(const CSolver* solver) { return sqrt(result); } +static void ApplyLineImplicitResidualSmoothing(CSolver* solver, CGeometry* geometry, unsigned short iMesh) { + if (iMesh == 0) return; + + const auto nPoint = geometry->GetnPointDomain(); + if (nPoint < 3) return; + + const unsigned short nVar = solver->GetnVar(); + const unsigned short nDim = geometry->GetnDim(); + const su2double damping = 0.25; + std::vector visited(nPoint, false); + unsigned long nLines = 0; + unsigned long totalLineLength = 0; + su2double totalLineResidual = 0.0; + + for (auto iSeed = 0ul; iSeed < nPoint; ++iSeed) { + if (visited[iSeed]) continue; + + std::vector line; + line.reserve(16); + line.push_back(iSeed); + visited[iSeed] = true; + + unsigned long current = iSeed; + for (int step = 0; step < 8; ++step) { + unsigned long next = std::numeric_limits::max(); + const auto* coordCurrent = geometry->nodes->GetCoord(current); + const unsigned short nNeigh = geometry->nodes->GetnPoint(current); + su2double bestScore = -1e30; + + for (unsigned short iNeigh = 0; iNeigh < nNeigh; ++iNeigh) { + const auto candidate = geometry->nodes->GetPoint(current, iNeigh); + if (candidate >= nPoint || candidate == current || visited[candidate]) continue; + + const auto* coordCandidate = geometry->nodes->GetCoord(candidate); + su2double score = 0.0; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) { + const su2double delta = fabs(coordCandidate[iDim] - coordCurrent[iDim]); + score += delta; + } + if (score > bestScore) { + bestScore = score; + next = candidate; + } + } + + if (next == std::numeric_limits::max()) break; + line.push_back(next); + visited[next] = true; + current = next; + } + + const auto nLine = static_cast(line.size()); + if (nLine < 2) continue; + + ++nLines; + totalLineLength += nLine; + + std::vector lineAverage(nVar, 0.0); + for (auto i = 0ul; i < nLine; ++i) { + const auto* residual = solver->GetNodes()->GetResidual_Old(line[i]); + if (residual == nullptr) continue; + for (unsigned short iVar = 0; iVar < nVar; ++iVar) { + lineAverage[iVar] += residual[iVar]; + totalLineResidual += fabs(residual[iVar]); + } + } + + for (unsigned short iVar = 0; iVar < nVar; ++iVar) { + lineAverage[iVar] /= static_cast(nLine); + } + + for (auto i = 0ul; i < nLine; ++i) { + const auto* oldResidual = solver->GetNodes()->GetResidual_Old(line[i]); + std::vector block(nVar, 0.0); + for (unsigned short iVar = 0; iVar < nVar; ++iVar) { + block[iVar] = oldResidual[iVar] + damping * (lineAverage[iVar] - oldResidual[iVar]); + } + solver->LinSysRes.SetBlock(line[i], block.data()); + } + } + + if (SU2_MPI::GetRank() == MASTER_NODE) { + const su2double avgLineLength = (nLines > 0) ? static_cast(totalLineLength) / static_cast(nLines) : 0.0; + const su2double avgLineResidual = (nLines > 0) ? totalLineResidual / static_cast(nLines) : 0.0; + cout << "[MG LINE] level=" << iMesh << " activated=" << (nLines > 0 ? "yes" : "no") + << " lines=" << nLines << " avgLen=" << avgLineLength + << " avgResidual=" << avgLineResidual << std::endl; + } +} + } // anonymous namespace void CMultiGridIntegration::adaptDampingFactors(CConfig* config, passivedouble crossCycleRatio) { @@ -248,49 +339,27 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, config[iZone]->GetnMGLevels()); } - /*--- Perform the Full Approximation Scheme multigrid ---*/ - - MultiGrid_Cycle(geometry, solver_container, numerics_container, config, - FinestMesh, RecursiveParam, RunTime_EqSystem, iZone, iInst); - - /*--- After a turb FAS V-cycle: recompute mu_t at the finest active level from the updated - * nu_hat/k/omega and restrict it to all coarser levels. The flow FAS on the NEXT outer - * iteration uses these mu_t values at every coarse level for the eddy-viscosity coupling. - * (Postprocessing was already called on FinestMesh inside the last PreSmoothing step of - * MultiGrid_Cycle; we call it once more to be safe after the V-cycle correction is applied.) ---*/ - if (RunTime_EqSystem == RUNTIME_TURB_SYS && - config[iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE) { - solver_container[iZone][iInst][FinestMesh][TURB_SOL]->Postprocessing( - geometry[iZone][iInst][FinestMesh], - solver_container[iZone][iInst][FinestMesh], - config[iZone], FinestMesh); - RestrictTurbEddyViscToCoarseLevels(geometry[iZone][iInst], - solver_container[iZone][iInst], - config[iZone], FinestMesh, - config[iZone]->GetnMGLevels()); - } - - /*--- Adapt coarse-grid CFL once per cycle using smoothing residuals gathered during the cycle. ---*/ + /*--- Rebuild coarse-grid CFL before the cycle so the currently active FMG + * level uses the intended CFL in this iteration. ---*/ const unsigned short nMGLevels = config[iZone]->GetnMGLevels(); BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - /*--- Use the current finest-grid CFL as the base for deterministic - * coarse-level scaling. Fall back to config scalar when local CFL - * adaptation is disabled. ---*/ + /*--- Use the level-0 flow CFL as the base reference and derive all coarse + * levels from it via MG_CFL_SCALING[i] = CFL(i+1)/CFL(i). Fall back to + * config scalar when local level-0 CFL is unavailable. ---*/ passivedouble cfl_base = SU2_TYPE::GetValue( - solver_container[iZone][iInst][FinestMesh][Solver_Position]->GetAvg_CFL_Local()); + solver_container[iZone][iInst][MESH_0][Solver_Position]->GetAvg_CFL_Local()); if (cfl_base < EPS) - cfl_base = SU2_TYPE::GetValue(config[iZone]->GetCFL(FinestMesh)); + cfl_base = SU2_TYPE::GetValue(config[iZone]->GetCFL(MESH_0)); const auto& cflScaling = config[iZone]->GetMGOptions().MG_CflScaling; passivedouble CFL_local = cfl_base; - for (unsigned short iMesh = FinestMesh; iMesh < nMGLevels; ++iMesh) { - const unsigned short lvl = iMesh + 1; - /*--- Use per-level scaling factor to increase coarse CFL (allows values > 1.0). - * Index into cflScaling is iMesh (0-based transition). ---*/ - const passivedouble scale = (iMesh < cflScaling.size()) - ? max(passivedouble{1e-6}, SU2_TYPE::GetValue(cflScaling[iMesh])) + for (unsigned short lvl = 1; lvl <= nMGLevels; ++lvl) { + /*--- Index into cflScaling is (lvl-1): transition lvl-1 -> lvl. ---*/ + const unsigned short iScale = lvl - 1; + const passivedouble scale = (iScale < cflScaling.size()) + ? max(passivedouble{1e-6}, SU2_TYPE::GetValue(cflScaling[iScale])) : passivedouble{0.25}; CFL_local *= scale; config[iZone]->SetCFL(lvl, CFL_local); @@ -298,17 +367,40 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, } END_SU2_OMP_SAFE_GLOBAL_ACCESS - /*--- Propagate the updated coarse-grid CFL to every coarse-grid point (all threads). ---*/ - for (unsigned short iMesh = FinestMesh; iMesh < nMGLevels; ++iMesh) { - const passivedouble CFL_coarse_new = SU2_TYPE::GetValue(config[iZone]->GetCFL(iMesh+1)); - CGeometry* geo_c = geometry[iZone][iInst][iMesh+1]; - CSolver* sol_c = solver_container[iZone][iInst][iMesh+1][Solver_Position]; + /*--- Propagate updated CFL to all coarse-grid points before the cycle. ---*/ + for (unsigned short iMesh = 1; iMesh <= nMGLevels; ++iMesh) { + const passivedouble CFL_coarse_new = SU2_TYPE::GetValue(config[iZone]->GetCFL(iMesh)); + CGeometry* geo_c = geometry[iZone][iInst][iMesh]; + CSolver* sol_c = solver_container[iZone][iInst][iMesh][Solver_Position]; + SU2_OMP_SAFE_GLOBAL_ACCESS(sol_c->SetCFL_Local_Stats(CFL_coarse_new);) SU2_OMP_FOR_STAT(roundUpDiv(geo_c->GetnPoint(), omp_get_num_threads())) for (auto iPoint = 0ul; iPoint < geo_c->GetnPoint(); iPoint++) sol_c->GetNodes()->SetLocalCFL(iPoint, CFL_coarse_new); END_SU2_OMP_FOR } + /*--- Perform the Full Approximation Scheme multigrid ---*/ + + MultiGrid_Cycle(geometry, solver_container, numerics_container, config, + FinestMesh, RecursiveParam, RunTime_EqSystem, iZone, iInst); + + /*--- After a turb FAS V-cycle: recompute mu_t at the finest active level from the updated + * nu_hat/k/omega and restrict it to all coarser levels. The flow FAS on the NEXT outer + * iteration uses these mu_t values at every coarse level for the eddy-viscosity coupling. + * (Postprocessing was already called on FinestMesh inside the last PreSmoothing step of + * MultiGrid_Cycle; we call it once more to be safe after the V-cycle correction is applied.) ---*/ + if (RunTime_EqSystem == RUNTIME_TURB_SYS && + config[iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE) { + solver_container[iZone][iInst][FinestMesh][TURB_SOL]->Postprocessing( + geometry[iZone][iInst][FinestMesh], + solver_container[iZone][iInst][FinestMesh], + config[iZone], FinestMesh); + RestrictTurbEddyViscToCoarseLevels(geometry[iZone][iInst], + solver_container[iZone][iInst], + config[iZone], FinestMesh, + config[iZone]->GetnMGLevels()); + } + /*--- Computes primitive variables and gradients in the finest mesh (useful for the next solver (turbulence) and output ---*/ solver_container[iZone][iInst][MESH_0][Solver_Position]->Preprocessing(geometry[iZone][iInst][MESH_0], @@ -413,8 +505,8 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, for (unsigned short i = 0; i <= nMGLevels; ++i) { std::ostringstream ss; ss << std::fixed << std::setprecision(4); - if (i == MESH_0) { - ss << SU2_TYPE::GetValue(solver_container[iZone][iInst][MESH_0][Solver_Position]->GetAvg_CFL_Local()); + if (i == FinestMesh) { + ss << SU2_TYPE::GetValue(solver_container[iZone][iInst][FinestMesh][Solver_Position]->GetAvg_CFL_Local()); } else { ss << SU2_TYPE::GetValue(config[iZone]->GetCFL(i)); } @@ -889,6 +981,10 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ } END_SU2_OMP_FOR + if (iMesh > 0 && RunTime_EqSystem == RUNTIME_FLOW_SYS) { + ApplyLineImplicitResidualSmoothing(solver, geometry, iMesh); + } + /*--- Restore original residuals (without average) at boundary points. ---*/ for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); iMarker++) { diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index bdfd71fb8fc..a1b67a4c8ef 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -139,11 +139,11 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe RUNTIME_RADIATION_SYS, val_iZone, val_iInst); } - /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. - During Full-MG warmup (FinestMesh > MESH_0), skip adaptation entirely until the finest - mesh is active. ---*/ + /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. + For Full-MG, allow adaptation during warmup as well so CFL can evolve on coarse active + levels before reaching the finest mesh. ---*/ SU2_OMP_PARALLEL - if (!disc_adj && config[val_iZone]->GetFinestMesh() == MESH_0) { + if (!disc_adj) { solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone]); solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->IdentifySolutionOutliers(config[val_iZone], InnerIter); @@ -176,7 +176,6 @@ void CFluidIteration::Update(COutput* output, CIntegration**** integration, CGeo if ((config[val_iZone]->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || (config[val_iZone]->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) { /*--- Update dual time solver on all mesh levels ---*/ - for (unsigned short iMesh = 0; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) { integration[val_iZone][val_iInst][FLOW_SOL]->SetDualTime_Solver(geometry[val_iZone][val_iInst][iMesh], solver[val_iZone][val_iInst][iMesh][FLOW_SOL], From 4f6def7db5163f962c6b4f59c4f1bd8c5fc4ccc7 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Tue, 4 Aug 2026 20:57:01 +0200 Subject: [PATCH 7/7] remove turbulence freezing --- Common/include/option_structure.hpp | 1 + Common/src/CConfig.cpp | 2 + .../integration/CMultiGridIntegration.hpp | 12 ++ .../integration/CSingleGridIntegration.hpp | 12 ++ SU2_CFD/include/solvers/CTurbSASolver.hpp | 19 +++ SU2_CFD/include/variables/CTurbVariable.hpp | 30 +++++ SU2_CFD/include/variables/CVariable.hpp | 12 ++ .../src/integration/CMultiGridIntegration.cpp | 117 ++++++++++++++++-- .../integration/CSingleGridIntegration.cpp | 49 ++++++++ SU2_CFD/src/iteration/CFluidIteration.cpp | 42 +++++++ SU2_CFD/src/solvers/CSolver.cpp | 4 +- SU2_CFD/src/solvers/CTurbSASolver.cpp | 86 +++++++++++++ SU2_CFD/src/variables/CTurbVariable.cpp | 2 + 13 files changed, 379 insertions(+), 9 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 14e09e9158f..df18786efd9 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1129,6 +1129,7 @@ struct CMGOptions { bool MG_Implicit_Lines_Isotropic{false}; /*!< \brief Use isotropic (vs anisotropic) agglomeration along implicit lines. */ unsigned long MG_Startup_Iter{100}; /*!< \brief Number of iterations on coarsest mesh during FMG startup phase. */ bool TurbMG{false}; /*!< \brief Run turbulence equations through a FAS MG V-cycle instead of single-grid. */ + bool MG_Turb_Freeze_Source{false}; /*!< \brief Freeze turbulence source terms on coarse multigrid levels. */ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 208767d6db5..adaf00bc671 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2075,6 +2075,8 @@ void CConfig::SetConfig_Options() { addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100); /*!\brief MG_TURB\n DESCRIPTION: Run turbulence equations through a FAS Multigrid V-cycle instead of single-grid. DEFAULT: NO \ingroup Config*/ addBoolOption("MG_TURB", MGOptions.TurbMG, false); + /*!\brief MG_TURB_FREEZE_SOURCE\n DESCRIPTION: Freeze turbulence source terms on coarse multigrid levels using values from the fine grid. Reduces stiffness on coarse grids. DEFAULT: NO \ingroup Config*/ + addBoolOption("MG_TURB_FREEZE_SOURCE", MGOptions.MG_Turb_Freeze_Source, false); /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i). If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/ addDoubleListOption("MG_CFL_SCALING", nMG_CflScaling_p, MG_CflScaling_p); diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp index 1b5c7476bf6..93a6d58a890 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -231,6 +231,18 @@ class CMultiGridIntegration final : public CIntegration { void SetRestricted_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + /*! + * \brief Restrict frozen turbulence source terms from fine grid to coarse grid. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] sol_fine - Pointer to the solution on the fine grid. + * \param[out] sol_coarse - Pointer to the solution on the coarse grid. + * \param[in] geo_fine - Geometrical definition of the fine grid. + * \param[in] geo_coarse - Geometrical definition of the coarse grid. + * \param[in] config - Definition of the particular problem. + */ + void SetRestricted_FrozenSource(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + /*! * \brief Initialize the adjoint solution using the primal problem. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/integration/CSingleGridIntegration.hpp b/SU2_CFD/include/integration/CSingleGridIntegration.hpp index 2c08096e969..2c527ad35e1 100644 --- a/SU2_CFD/include/integration/CSingleGridIntegration.hpp +++ b/SU2_CFD/include/integration/CSingleGridIntegration.hpp @@ -63,6 +63,18 @@ class CSingleGridIntegration final : public CIntegration { void SetRestricted_EddyVisc(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + /*! + * \brief Restrict frozen turbulence source terms from fine grid to coarse grid. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] sol_fine - Pointer to the solution on the fine grid. + * \param[out] sol_coarse - Pointer to the solution on the coarse grid. + * \param[in] geo_fine - Geometrical definition of the fine grid. + * \param[in] geo_coarse - Geometrical definition of the coarse grid. + * \param[in] config - Definition of the particular problem. + */ + void SetRestricted_FrozenSource(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + public: /*! * \brief Constructor of the class. diff --git a/SU2_CFD/include/solvers/CTurbSASolver.hpp b/SU2_CFD/include/solvers/CTurbSASolver.hpp index fa2cc654019..1b22b4985c9 100644 --- a/SU2_CFD/include/solvers/CTurbSASolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSASolver.hpp @@ -43,6 +43,25 @@ class CTurbSASolver final : public CTurbSolver { su2double nu_tilde_Engine[4] = {0.0}; su2double nu_tilde_ActDisk[4] = {0.0}; + /*! + * \brief Override SetTime_Step to include source term stiffness in the turbulence time step. + * \details The base CScalarSolver implementation scales the flow time step by the CFL ratio, + * which only captures convective physics. For the SA model the source term destruction + * dS/d(nu_tilde) ~ c_w*nu_tilde/d^2 dominates near walls (d->0), adding large negative + * diagonal contributions. The frozen source Jacobian (cached each fine-grid iteration) + * captures this stiffness. Limiting dt so that the diagonal is source-dominated rather + * than time-dominated prevents the implicit update from being O(R*dt/V->0), which + * produces r=1.000 (no residual reduction) in the multigrid pre-smoother. + * + * \param[in] geometry - Geometrical definition. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] Iteration - External iteration number. + */ + void SetTime_Step(CGeometry* geometry, CSolver** solver_container, CConfig* config, + unsigned short iMesh, unsigned long Iteration) override; + /*! * \brief A virtual member. * \param[in] solver - Solver container diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index f10ebe5f1a1..74a543b4d2e 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -38,6 +38,8 @@ class CTurbVariable : public CScalarVariable { protected: VectorType muT; /*!< \brief Eddy viscosity. */ + VectorType frozen_source; /*!< \brief Frozen source term from fine grid for multigrid coarse levels. */ + VectorType frozen_source_jacobian; /*!< \brief Frozen source Jacobian diagonal for implicit coupling on coarse grids. */ public: static constexpr size_t MAXNVAR = 4; @@ -106,6 +108,34 @@ class CTurbVariable : public CScalarVariable { * \param[in] val_DC_kw - diffusion coefficient value */ + /*! + * \brief Get the frozen source term for multigrid coarse levels. + * \param[in] iPoint - Point index. + * \return Frozen source term. + */ + inline su2double GetFrozenSource(unsigned long iPoint) const { return frozen_source(iPoint); } + + /*! + * \brief Set the frozen source term for multigrid coarse levels. + * \param[in] iPoint - Point index. + * \param[in] val_source - Frozen source term value. + */ + inline void SetFrozenSource(unsigned long iPoint, su2double val_source) { frozen_source(iPoint) = val_source; } + + /*! + * \brief Get the frozen source Jacobian diagonal for multigrid coarse levels. + * \param[in] iPoint - Point index. + * \return Frozen source Jacobian diagonal. + */ + inline su2double GetFrozenSourceJacobian(unsigned long iPoint) const { return frozen_source_jacobian(iPoint); } + + /*! + * \brief Set the frozen source Jacobian diagonal for multigrid coarse levels. + * \param[in] iPoint - Point index. + * \param[in] val_jac - Frozen source Jacobian diagonal value. + */ + inline void SetFrozenSourceJacobian(unsigned long iPoint, su2double val_jac) { frozen_source_jacobian(iPoint) = val_jac; } + /*! * \brief Register eddy viscosity (muT) as Input or Output of an AD recording. * \param[in] input - Boolean whether In- or Output should be registered. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 528adc139db..12cde687fd6 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1778,6 +1778,18 @@ class CVariable { */ inline virtual su2double GetTurbIndex(unsigned long iPoint) const {return 0.0;} + /*! + * \brief Get the frozen turbulence source term density for multigrid. + * \return Frozen source density (source per unit volume). + */ + inline virtual su2double GetFrozenSource(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief Set the frozen turbulence source term density for multigrid. + * \param[in] val_source - Frozen source density value. + */ + inline virtual void SetFrozenSource(unsigned long iPoint, su2double val_source) {} + /*! * \brief A virtual member. * \param[in] iVar - Index of the variable. diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 2429a0e33d1..53a686e0734 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -28,6 +28,7 @@ #include "../../include/integration/CMultiGridIntegration.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" +#include "../../include/variables/CTurbVariable.hpp" #include #include #include @@ -594,10 +595,21 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, /*--- Assemble the fine-grid defect term that will be restricted to the coarse-grid FAS problem. ---*/ SetResidual_Term(geometry_fine, solver_fine); + /*--- Communicate frozen sources on finest level to halo cells before restriction aggregates from fine children. ---*/ + if (iMesh == MESH_0 && RunTime_EqSystem == RUNTIME_TURB_SYS && config->GetMGOptions().MG_Turb_Freeze_Source) { + solver_fine->InitiateComms(geometry_fine, config, MPI_QUANTITIES::SOLUTION_EDDY); + solver_fine->CompleteComms(geometry_fine, config, MPI_QUANTITIES::SOLUTION_EDDY); + } + /*--- Restrict the fine-grid state to the coarse grid and initialize the coarse-grid state. ---*/ SetRestricted_Solution(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); + /*--- Restrict frozen source terms for turbulence multigrid if feature is enabled. ---*/ + if (RunTime_EqSystem == RUNTIME_TURB_SYS && config->GetMGOptions().MG_Turb_Freeze_Source) { + SetRestricted_FrozenSource(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); + } + solver_coarse->Preprocessing(geometry_coarse, solver_container_coarse, config, iMesh+1, NO_RK_ITER, RunTime_EqSystem, false); /*--- For turbulence: ensure flow primitives (density, laminar viscosity) are updated on the @@ -887,7 +899,7 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS /*--- Enforce Euler wall BC on corrections by projecting to tangent plane ---*/ sol_coarse->MultigridProjectEulerWall(geo_coarse, config, true); - /*--- Remove any contributions from no-slip walls. ---*/ + /*--- Remove any contributions from no-slip walls (Dirichlet BC enforcement). ---*/ for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetViscous_Wall(iMarker)) { @@ -896,11 +908,19 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { auto Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode(); - /*--- For dirichlet boundary conditions, set the correction to zero. - Note that Solution_Old stores the correction not the actual value ---*/ + /*--- For Dirichlet boundary conditions, set the correction to zero. + * Note that Solution_Old stores the correction, not the actual value. ---*/ - su2double zero[3] = {0.0}; - sol_coarse->GetNodes()->SetVelocity_Old(Point_Coarse, zero); + if (RunTime_EqSystem == RUNTIME_TURB_SYS) { + /*--- Turbulence: explicitly zero all variables (scalar equations). ---*/ + for (auto iVar = 0u; iVar < sol_coarse->GetnVar(); iVar++) { + sol_coarse->GetNodes()->SetSolution_Old(Point_Coarse, iVar, 0.0); + } + } else { + /*--- Flow: zero velocity components only (momentum equations). ---*/ + su2double zero[3] = {0.0}; + sol_coarse->GetNodes()->SetVelocity_Old(Point_Coarse, zero); + } } END_SU2_OMP_FOR @@ -975,7 +995,13 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ for (auto iVar = 0u; iVar < nVar; iVar++) { su2double smoothed = (Residual_Old[iVar] + val_smooth_coeff*Residual_Sum[iVar])*factor; - if (use_conservative_damping) smoothed *= turbulence_base_damping; + /*--- FIX 2: Do NOT apply turbulence_base_damping (0.50) to turbulence corrections. + * use_conservative_damping is true for nVar<=2 which includes SA (nVar=1), but the + * 0.50 factor was designed for compressible flow density-velocity-pressure coupling + * near walls. Applying it to a scalar equation halves all corrections unnecessarily. + * Combined with SetProlongated_Correction's damping (~0.375 at Level 1) and the config + * damping factor (0.9), only 17% of the coarse correction was reaching the fine grid. ---*/ + if (use_conservative_damping && RunTime_EqSystem == RUNTIME_FLOW_SYS) smoothed *= turbulence_base_damping; solver->LinSysRes(iPoint,iVar) = smoothed; } } @@ -1021,7 +1047,12 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet SU2_ZONE_SCOPED const unsigned short nVar = sol_fine->GetnVar(); - const bool use_conservative_damping = (nVar <= 2); + /*--- Conservative damping prevents oscillations from density-velocity-pressure coupling + * in compressible flow (nVar > 2 equations). Scalar/turbulence equations (nVar <= 2) + * are single-variable and have no such coupling; applying flow-level damping (base ~0.375) + * reduces turbulence corrections to only 34% of their computed value, preventing the + * coarse-grid MG cycle from being effective. ---*/ + const bool use_conservative_damping = (nVar > 2); const su2double levelScale = GetMGLevelCorrectionScale(iMesh); const su2double base_damping = use_conservative_damping ? max(su2double{0.15}, 0.50 * levelScale) : 1.0; const su2double wall_damping = use_conservative_damping ? max(su2double{0.10}, 0.25 * levelScale) : 1.0; @@ -1201,7 +1232,18 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar SU2_OMP_FOR_STAT(32) for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { auto Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode(); - sol_coarse->GetNodes()->SetVel_ResTruncError_Zero(Point_Coarse); + /*--- FIX 1: For turbulence (Dirichlet BC nu_tilde=0 at walls), zero the FULL truncation error. + * SetVel_ResTruncError_Zero is a no-op for CTurbVariable (only zeros velocity components + * in flow variables). Without this, the FAS forcing term injects spurious residuals at + * wall nodes where the coarse grid should simply enforce nu_tilde=0. For flow solvers, + * only the velocity components need zeroing (pressure and energy BCs are flux-based). ---*/ + if (sol_coarse->GetnVar() <= 2) { + /*--- Turbulence/scalar: zero ALL truncation error components at the Dirichlet wall. ---*/ + sol_coarse->GetNodes()->SetRes_TruncErrorZero(Point_Coarse); + } else { + /*--- Flow: zero only velocity components (pressure/energy BCs are flux-based). ---*/ + sol_coarse->GetNodes()->SetVel_ResTruncError_Zero(Point_Coarse); + } } END_SU2_OMP_FOR } @@ -1295,6 +1337,17 @@ void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSyst sol_coarse->GetNodes()->SetVelSolutionDVector(Point_Coarse); } + if (Solver_Position == TURB_SOL) { + /*--- CRITICAL FIX: Enforce Dirichlet BC for turbulence at walls. + * SA model requires nu_tilde = 0 at smooth walls. After restriction, + * wall values are averaged from fine grid children, violating the BC. + * This causes coarse grid solves with incorrect boundary conditions, + * producing corrupted corrections that stall convergence. ---*/ + for (auto iVar = 0u; iVar < sol_coarse->GetnVar(); iVar++) { + sol_coarse->GetNodes()->SetSolution(Point_Coarse, iVar, 0.0); + } + } + } END_SU2_OMP_FOR } @@ -1310,6 +1363,54 @@ void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSyst } +void CMultiGridIntegration::SetRestricted_FrozenSource(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { + SU2_ZONE_SCOPED + + /*--- Only applicable for turbulence equations ---*/ + if (RunTime_EqSystem != RUNTIME_TURB_SYS) return; + + auto* turbNodes_fine = su2staticcast_p(sol_fine->GetNodes()); + auto* turbNodes_coarse = su2staticcast_p(sol_coarse->GetNodes()); + + /*--- Compute coarse frozen source density from fine grid using volume-weighted averaging. + * This is analogous to eddy viscosity restriction - both are intensive properties. ---*/ + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + + su2double Volume_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); + su2double SourceDensity_Coarse = 0.0; + su2double SourceJacobian_Coarse = 0.0; + + /*--- Volume-weighted average of source density and Jacobian from fine grid children ---*/ + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + su2double Volume_Child = geo_fine->nodes->GetVolume(Point_Fine); + su2double weight = Volume_Child / Volume_Parent; + + su2double SourceDensity_Fine = turbNodes_fine->GetFrozenSource(Point_Fine); + SourceDensity_Coarse += SourceDensity_Fine * weight; + + /*--- Restrict Jacobian (intensive property, volume-weighted averaging) ---*/ + su2double SourceJac_Fine = turbNodes_fine->GetFrozenSourceJacobian(Point_Fine); + SourceJacobian_Coarse += SourceJac_Fine * weight; + } + + turbNodes_coarse->SetFrozenSource(Point_Coarse, SourceDensity_Coarse); + turbNodes_coarse->SetFrozenSourceJacobian(Point_Coarse, SourceJacobian_Coarse); + } + END_SU2_OMP_FOR + + /*--- MPI communication of restricted frozen source to halo cells. + * Although frozen_source is not part of the solution vector, halo cells need correct values + * for proper parallel execution, especially when coarse points near partition boundaries + * aggregate from fine grid children. We reuse SOLUTION_EDDY communication infrastructure. ---*/ + + sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_EDDY); + sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_EDDY); + +} + void CMultiGridIntegration::SetRestricted_Gradient(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { SU2_ZONE_SCOPED diff --git a/SU2_CFD/src/integration/CSingleGridIntegration.cpp b/SU2_CFD/src/integration/CSingleGridIntegration.cpp index fab97842048..736d00ef2bb 100644 --- a/SU2_CFD/src/integration/CSingleGridIntegration.cpp +++ b/SU2_CFD/src/integration/CSingleGridIntegration.cpp @@ -26,6 +26,7 @@ */ #include "../../include/integration/CSingleGridIntegration.hpp" +#include "../../include/variables/CTurbVariable.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" @@ -100,6 +101,16 @@ void CSingleGridIntegration::SingleGrid_Iteration(CGeometry ****geometry, CSolve geometry[iZone][iInst][iMesh], geometry[iZone][iInst][iMesh+1], config[iZone]); + + /*--- Restrict frozen source terms if feature is enabled ---*/ + if (config[iZone]->GetMGOptions().MG_Turb_Freeze_Source) { + SetRestricted_FrozenSource(RunTime_EqSystem, + solver_container[iZone][iInst][iMesh][Solver_Position], + solver_container[iZone][iInst][iMesh+1][Solver_Position], + geometry[iZone][iInst][iMesh], + geometry[iZone][iInst][iMesh+1], + config[iZone]); + } } } @@ -168,3 +179,41 @@ void CSingleGridIntegration::SetRestricted_EddyVisc(unsigned short RunTime_EqSys sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_EDDY); } + +void CSingleGridIntegration::SetRestricted_FrozenSource(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { + SU2_ZONE_SCOPED + + /*--- Cast to turbulence variable type to access frozen source methods ---*/ + auto* turbNodes_fine = su2staticcast_p(sol_fine->GetNodes()); + auto* turbNodes_coarse = su2staticcast_p(sol_coarse->GetNodes()); + + /*--- Compute coarse frozen source density from fine grid using volume-weighted averaging. + * This is analogous to eddy viscosity restriction - both are intensive properties. ---*/ + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + + su2double Volume_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); + su2double SourceDensity_Coarse = 0.0; + + /*--- Volume-weighted average of source density from fine grid children ---*/ + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + su2double Volume_Child = geo_fine->nodes->GetVolume(Point_Fine); + su2double SourceDensity_Fine = turbNodes_fine->GetFrozenSource(Point_Fine); + SourceDensity_Coarse += SourceDensity_Fine * Volume_Child / Volume_Parent; + } + + turbNodes_coarse->SetFrozenSource(Point_Coarse, SourceDensity_Coarse); + } + END_SU2_OMP_FOR + + /*--- MPI communication of restricted frozen source to halo cells. + * Although frozen_source is not part of the solution vector, halo cells need correct values + * for proper parallel execution, especially when coarse points near partition boundaries + * aggregate from fine grid children. We reuse SOLUTION_EDDY communication infrastructure. ---*/ + + sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_EDDY); + sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_EDDY); + +} diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index a1b67a4c8ef..b2fa08845bb 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -100,8 +100,50 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_TURB_SYS); if (config[val_iZone]->GetMGOptions().TurbMG) { + /*--- Before turbulence multigrid: restrict updated FLOW solution from MESH_0 to all coarse grids. + * After flow MG completes, coarse grids have stale flow data (only MESH_0 is updated via prolongation). + * Turbulence MG needs current flow primitives (density, laminar viscosity) on coarse grids for source terms. + * Without this restriction, coarse-grid FLOW Preprocessing uses stale conservative variables, producing + * incorrect primitives that corrupt eddy viscosity computation and cause oscillations. ---*/ + for (auto iMesh = 1u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) { + CSolver::MultigridRestriction(*geometry[val_iZone][val_iInst][iMesh - 1], + solver[val_iZone][val_iInst][iMesh - 1][FLOW_SOL]->GetNodes()->GetSolution(), + *geometry[val_iZone][val_iInst][iMesh], + solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->GetNodes()->GetSolution()); + /*--- Synchronize halo cells across MPI ranks before preprocessing (critical for parallel runs) ---*/ + solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->InitiateComms(geometry[val_iZone][val_iInst][iMesh], + config[val_iZone], MPI_QUANTITIES::SOLUTION); + solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->CompleteComms(geometry[val_iZone][val_iInst][iMesh], + config[val_iZone], MPI_QUANTITIES::SOLUTION); + /*--- Update flow primitives on coarse grid from newly restricted conservative variables ---*/ + solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->Preprocessing(geometry[val_iZone][val_iInst][iMesh], + solver[val_iZone][val_iInst][iMesh], + config[val_iZone], iMesh, NO_RK_ITER, + RUNTIME_FLOW_SYS, false); + } + integration[val_iZone][val_iInst][TURB_SOL]->MultiGrid_Iteration(geometry, solver, numerics, config, RUNTIME_TURB_SYS, val_iZone, val_iInst); + + /*--- After turbulence multigrid: restrict updated TURBULENCE solution from MESH_0 to all coarse grids. + * This ensures coarse grids have current turbulence values for the next iteration's FLOW MG. + * Without this, coarse-grid flow equations use stale eddy viscosity, creating a feedback loop + * that eventually causes oscillations after accumulated error builds up over many iterations. ---*/ + for (auto iMesh = 1u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) { + CSolver::MultigridRestriction(*geometry[val_iZone][val_iInst][iMesh - 1], + solver[val_iZone][val_iInst][iMesh - 1][TURB_SOL]->GetNodes()->GetSolution(), + *geometry[val_iZone][val_iInst][iMesh], + solver[val_iZone][val_iInst][iMesh][TURB_SOL]->GetNodes()->GetSolution()); + /*--- Synchronize halo cells across MPI ranks before postprocessing (critical for parallel runs) ---*/ + solver[val_iZone][val_iInst][iMesh][TURB_SOL]->InitiateComms(geometry[val_iZone][val_iInst][iMesh], + config[val_iZone], MPI_QUANTITIES::SOLUTION); + solver[val_iZone][val_iInst][iMesh][TURB_SOL]->CompleteComms(geometry[val_iZone][val_iInst][iMesh], + config[val_iZone], MPI_QUANTITIES::SOLUTION); + /*--- Update turbulence postprocessing (eddy viscosity) on coarse grid ---*/ + solver[val_iZone][val_iInst][iMesh][TURB_SOL]->Postprocessing(geometry[val_iZone][val_iInst][iMesh], + solver[val_iZone][val_iInst][iMesh], + config[val_iZone], iMesh); + } } else { integration[val_iZone][val_iInst][TURB_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_TURB_SYS, val_iZone, val_iInst); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index e014c794a52..1e073bf16ef 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1365,7 +1365,7 @@ void CSolver::GetCommCountAndType(const CConfig* config, MPI_TYPE = COMM_TYPE::DOUBLE; break; case MPI_QUANTITIES::SOLUTION_EDDY: - COUNT_PER_POINT = nVar+1; + COUNT_PER_POINT = nVar+2; MPI_TYPE = COMM_TYPE::DOUBLE; break; case MPI_QUANTITIES::STOCH_SOURCE_LANG: @@ -1503,6 +1503,7 @@ void CSolver::InitiateComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar); bufDSend[buf_offset+nVar] = base_nodes->GetmuT(iPoint); + bufDSend[buf_offset+nVar+1] = base_nodes->GetFrozenSource(iPoint); break; case MPI_QUANTITIES::STOCH_SOURCE_LANG: for (iDim = 0; iDim < nDim; iDim++) @@ -1659,6 +1660,7 @@ void CSolver::CompleteComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]); base_nodes->SetmuT(iPoint,bufDRecv[buf_offset+nVar]); + base_nodes->SetFrozenSource(iPoint, bufDRecv[buf_offset+nVar+1]); break; case MPI_QUANTITIES::STOCH_SOURCE_LANG: for (iDim = 0; iDim < nDim; iDim++) diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 714d1972695..2c3a4ca3754 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -199,6 +199,24 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor } +void CTurbSASolver::SetTime_Step(CGeometry* geometry, CSolver** solver_container, CConfig* config, + unsigned short iMesh, unsigned long Iteration) { + SU2_ZONE_SCOPED + + const auto flowNodes = solver_container[FLOW_SOL]->GetNodes(); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + /*--- Scale flow time step by the turbulence CFL ratio (convective physics). + * The SA source Jacobian near walls is negative (destruction-dominated), which means + * SubtractBlock2Diag ADDS a positive contribution to the implicit diagonal — the source + * term improves conditioning near walls. No additional time step limiting is needed. ---*/ + const su2double dt = nodes->GetLocalCFL(iPoint) / flowNodes->GetLocalCFL(iPoint) * flowNodes->GetDelta_Time(iPoint); + nodes->SetDelta_Time(iPoint, dt); + } + END_SU2_OMP_FOR +} + void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { SU2_ZONE_SCOPED @@ -362,6 +380,12 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai const bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE); const bool transition_BC = config->GetSAParsedOptions().bc; + /*--- Check if we should use frozen source terms on coarse grids. ---*/ + const bool use_frozen_source = (iMesh > 0) && config->GetMGOptions().MG_Turb_Freeze_Source; + + /*--- Diagnostic counters for frozen vs computed source terms. ---*/ + unsigned long n_frozen = 0, n_computed = 0; + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); /*--- Pick one numerics object per thread. ---*/ @@ -374,6 +398,32 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai SU2_OMP_FOR_DYN(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + /*--- Check if we should use frozen source or compute it. ---*/ + if (use_frozen_source) { + /*--- Use frozen source density from fine grid (cached during restriction). ---*/ + su2double SourceDensity = nodes->GetFrozenSource(iPoint); + su2double Volume = geometry->nodes->GetVolume(iPoint); + + /*--- Add source contribution: LinSysRes += Source * Volume. ---*/ + LinSysRes(iPoint, 0) += SourceDensity * Volume; + + /*--- Add frozen source Jacobian to maintain implicit coupling. + * This is critical for convergence: without the Jacobian, the coarse grid + * linear system produces identical corrections every iteration, causing stalling. + * The Jacobian was cached on the fine grid and restricted here. ---*/ + if (implicit) { + su2double SourceJac = nodes->GetFrozenSourceJacobian(iPoint); + su2double jac_block[1][1]; + jac_block[0][0] = SourceJac * Volume; + Jacobian.SubtractBlock2Diag(iPoint, jac_block); + } + + SU2_OMP_ATOMIC + n_frozen++; + + } else { + /*--- Compute source term normally (fine grid or freezing disabled). ---*/ + /*--- Conservative variables w/o reconstruction ---*/ numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), nullptr); @@ -463,6 +513,27 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai if (implicit) Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); + /*--- Cache source density on fine grid when MG_TURB_FREEZE_SOURCE is enabled. + * The cached values will be restricted to coarse grids and used there instead + * of recomputing sources with stale flow data. ---*/ + if ((iMesh == 0) && config->GetMGOptions().MG_Turb_Freeze_Source) { + /*--- Store source density S (intensive property) for volume-weighted restriction. + * residual = -S*V from numerics, so S = -residual/V. + * This is analogous to storing eddy viscosity (also an intensive property). ---*/ + su2double Volume = geometry->nodes->GetVolume(iPoint); + nodes->SetFrozenSource(iPoint, -residual[0] / Volume); + + /*--- Cache source Jacobian diagonal for implicit coupling on coarse grids. + * Store as jacobian per unit volume (intensive property) for restriction. + * This preserves implicit coupling when using frozen sources. ---*/ + nodes->SetFrozenSourceJacobian(iPoint, residual.jacobian_i[0][0] / Volume); + } + + SU2_OMP_ATOMIC + n_computed++; + + } /* end of if-else use_frozen_source */ + } END_SU2_OMP_FOR @@ -490,6 +561,21 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai CustomSourceResidual(geometry, solver_container, numerics_container, config, iMesh); } + /*--- Diagnostic output for frozen source terms. ---*/ + if (config->GetMGOptions().MG_Turb_Freeze_Source) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + unsigned long n_frozen_global = 0, n_computed_global = 0; + SU2_MPI::Allreduce(&n_frozen, &n_frozen_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&n_computed, &n_computed_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + + if (SU2_MPI::GetRank() == MASTER_NODE && (config->GetInnerIter() % 100 == 0 || config->GetInnerIter() < 5)) { + cout << "[TURB_FREEZE] Mesh " << iMesh << ": Frozen=" << n_frozen_global + << " Computed=" << n_computed_global << endl; + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + } void CTurbSASolver::Source_Template(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp index 62b87f9363d..ede814b2d62 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -34,6 +34,8 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned turb_index.resize(nPoint) = su2double(1.0); intermittency.resize(nPoint) = su2double(1.0); + frozen_source.resize(nPoint) = su2double(0.0); + frozen_source_jacobian.resize(nPoint) = su2double(0.0); /*--- Allocate residual structures for multigrid (required for turbulence MG). ---*/ Res_TruncError.resize(nPoint, nVar) = su2double(0.0);