From 9a2d329e1b2ab456a061a15f49acf62155d46e46 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 29 Jul 2026 19:36:28 +0200 Subject: [PATCH 01/23] 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 02/23] 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 03/23] 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 04/23] 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 05/23] 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 3e0d528e4d2e0dd1896b400a33a203dd4afddd21 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 4 Aug 2026 21:28:47 +0200 Subject: [PATCH 06/23] Apply suggestions from code review remove multigrid turbulence Co-authored-by: Nijso --- Common/include/option_structure.hpp | 1 - Common/src/CConfig.cpp | 2 - .../integration/CMultiGridIntegration.hpp | 46 ----- .../src/integration/CMultiGridIntegration.cpp | 163 +----------------- SU2_CFD/src/iteration/CFluidIteration.cpp | 10 +- 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, 7 insertions(+), 237 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 14e09e9158f..e277e514332 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1128,7 +1128,6 @@ 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 208767d6db5..2d68a229b27 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2073,8 +2073,6 @@ 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 1b5c7476bf6..0672d1ae5ff 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -81,18 +81,6 @@ 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. @@ -192,20 +180,6 @@ 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. @@ -274,26 +248,6 @@ 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 12997293c87..0f9424191c6 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -184,46 +184,6 @@ 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();) } @@ -231,45 +191,11 @@ 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 @@ -469,8 +395,7 @@ 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); - /*--- 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. ---*/ +/*--- 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 ---*/ if ( iMesh < config->GetnMGLevels() ) { @@ -508,23 +433,10 @@ 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); - /*--- 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); - } + SetForcing_Term(solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, iMesh+1); /*--- Restore the time integration settings. ---*/ @@ -546,12 +458,7 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, iMesh+1, nextRecurseParam, RunTime_EqSystem, iZone, iInst); } - /*--- 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); - } + 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); @@ -912,11 +819,6 @@ 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"; - } } } @@ -1131,32 +1033,6 @@ 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 @@ -1350,36 +1226,3 @@ 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 bdfd71fb8fc..c3e95355858 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/multi-grid integrations. ---*/ + /*--- If the flow integration is not fully coupled, run the various single-grid integrations. ---*/ if (config[val_iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE && !frozen_visc) { @@ -95,14 +95,8 @@ 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); - - 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); } @@ -110,7 +104,6 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe 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); @@ -127,7 +120,6 @@ 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 00e7064d31c..ee798c1384b 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 = config->GetMGOptions().TurbMG ? INTEGRATION_TYPE::MULTIGRID : INTEGRATION_TYPE::SINGLEGRID; + metaData.integrationType = 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 714d1972695..b58afe4c9b3 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 || config->GetMGOptions().TurbMG) { + if (iMesh == MESH_0 || config->GetMGCycle() == MG_CYCLE::FULL) { /*--- 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 97b4ee1d78c..579b5f3cf72 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 || config->GetMGOptions().TurbMG) { + if (iMesh == MESH_0 || config->GetMGCycle() == MG_CYCLE::FULL) { /*--- 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 62b87f9363d..139223055a4 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -35,18 +35,6 @@ 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 2c7d882c11f..45a822755c9 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1732,10 +1732,6 @@ 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 4922e7b675398d8b1d93804dda3c1b554d5ebc26 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 4 Aug 2026 21:44:49 +0200 Subject: [PATCH 07/23] Apply suggestions from code review Co-authored-by: Nijso --- .../src/integration/CMultiGridIntegration.cpp | 80 ++----------------- SU2_CFD/src/iteration/CFluidIteration.cpp | 4 +- 2 files changed, 7 insertions(+), 77 deletions(-) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 0f9424191c6..74d2b2e75df 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -395,7 +395,7 @@ 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 ---*/ + /*--- 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 ---*/ if ( iMesh < config->GetnMGLevels() ) { @@ -418,16 +418,15 @@ 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 the fine-grid defect term is assembled. ---*/ + /*--- LinSysRes = R(u_N) here, before tau is added by SetResidual_Term. ---*/ 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); - /*--- Restrict the fine-grid state to the coarse grid and initialize the coarse-grid state. ---*/ + /*--- Compute $r_(k+1) = F_(k+1)(I^(k+1)_k u_k)$ ---*/ SetRestricted_Solution(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); @@ -435,6 +434,8 @@ 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); @@ -883,77 +884,6 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet } 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 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 - << "; 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; - 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 L" << iMesh << "] damp(wall/interior)= " << wall_damping << "/" << base_damping - << " raw max(wall/interior)= " << maxWall << "/" << maxInter - << " applied max(wall/interior)= " << maxWallApply << "/" << maxInterApply - << " levelScale=" << levelScale << "\n"; - } - } - END_SU2_OMP_SAFE_GLOBAL_ACCESS - } - /*--- MPI the new interpolated solution ---*/ sol_fine->InitiateComms(geo_fine, config, MPI_QUANTITIES::SOLUTION); diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index c3e95355858..89fee917a5c 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 grid integrations. ---*/ if (config[val_iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE && !frozen_visc) { @@ -97,7 +97,7 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe /*--- Solve the turbulence model ---*/ config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_TURB_SYS); - integration[val_iZone][val_iInst][TURB_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, + integration[val_iZone][val_iInst][TURB_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_TURB_SYS, val_iZone, val_iInst); } } From 5d806798b303aa07d99fa2dabe45e62f72ded981 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Tue, 4 Aug 2026 21:47:10 +0200 Subject: [PATCH 08/23] precommit --- Common/src/geometry/CMultiGridGeometry.cpp | 102 ++++++++++++--------- SU2_CFD/src/iteration/CFluidIteration.cpp | 1 - 2 files changed, 61 insertions(+), 42 deletions(-) diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index c9192e17afd..feeb1e2c60a 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -437,11 +437,16 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un 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 == 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++; @@ -459,9 +464,8 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un 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; + 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; } @@ -546,11 +550,16 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un 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 == 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++; @@ -559,9 +568,8 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un 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; + 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; } @@ -575,20 +583,25 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un 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++; + 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; + 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; + cout << " WARNING: Detected unexpected CV child counts (3-child=" << nCVs_3child << ", 4-child=" << nCVs_4child + << " in ANISO mode)" << endl; } } @@ -1378,9 +1391,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 ---*/ + 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 @@ -1484,7 +1497,8 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, 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; + 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; @@ -1597,8 +1611,8 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, 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 << "]" + 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; } @@ -1670,17 +1684,20 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, 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)); + 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; } + 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; + 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) ---*/ @@ -1697,12 +1714,15 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, 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; } + 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] << ")" + 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; } diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index 89fee917a5c..264effbdf7a 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -99,7 +99,6 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe 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]->GetKind_Species_Model() != SPECIES_MODEL::NONE) { From 7a0db1007efba030850b98a9978c9fcf2158b922 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 4 Aug 2026 21:50:27 +0200 Subject: [PATCH 09/23] Apply suggestions from code review Co-authored-by: Nijso --- SU2_CFD/src/integration/CMultiGridIntegration.cpp | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 74d2b2e75df..fae22207943 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -748,8 +748,6 @@ 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++) { @@ -789,11 +787,8 @@ 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++) { - 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; - } + for (auto iVar = 0u; iVar < nVar; iVar++) + solver->LinSysRes(iPoint,iVar) = (Residual_Old[iVar] + val_smooth_coeff*Residual_Sum[iVar])*factor; } END_SU2_OMP_FOR @@ -828,10 +823,6 @@ 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 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(); From aca34afd1a21a986112adc7a1952b67f6f570643 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 4 Aug 2026 21:55:35 +0200 Subject: [PATCH 10/23] Apply suggestions from code review Co-authored-by: Nijso --- .../src/integration/CMultiGridIntegration.cpp | 37 +------------------ 1 file changed, 2 insertions(+), 35 deletions(-) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index fae22207943..a14ad116f03 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -827,49 +827,16 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet /*--- 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]) { - 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; - constexpr su2double maxAllowedRatio = 1.25; - if (residualMag > 1e-30 && correctionMag > 1e-30) { - const su2double ratio = correctionMag / residualMag; - if (ratio > maxAllowedRatio) { - correctionScale = maxAllowedRatio / ratio; - } - } + if (Residual_Fine[iVar] != Residual_Fine[iVar]) +Residual_Fine[iVar] = 0.0; - 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; } } From 5648a3c54af211e7c435feb48fe6370a04daa9b0 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 09:26:10 +0200 Subject: [PATCH 11/23] fix cfl for full multigrid --- .../src/integration/CMultiGridIntegration.cpp | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index a14ad116f03..b873605108e 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -185,6 +185,23 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, config[iZone]); SU2_OMP_SAFE_GLOBAL_ACCESS(config[iZone]->SubtractFinestMesh();) + + auto nMG = config[iZone]->GetnMGLevels(); + auto cfl_base = SU2_TYPE::GetValue(config[iZone]->GetCFL(FinestMesh)); + auto cflScaling = config[iZone]->GetMGOptions().MG_CflScaling; + cout << "level =" << nMG - FinestMesh << endl; + cout <<"finestmesh = " << FinestMesh << endl; + cout << "CFL = " << cfl_base << ", scaling factor = " << cflScaling[FinestMesh-1] << endl; + // now scale the cfl using the scaling factors by multiplying from the first factor up to the factor for the current level + for (unsigned short iMesh = 0; iMesh < FinestMesh-1; ++iMesh) { + cout << "scale = " << cflScaling[iMesh] << endl; + cfl_base *= cflScaling[iMesh]; + } + cout << "new CFL = " << cfl_base << endl; + // set the new base cfl for MESH_0 to the scaled value: + SU2_OMP_SAFE_GLOBAL_ACCESS(config[iZone]->SetCFL(MESH_0, cfl_base);) + + } /*--- Set the current finest grid (full multigrid strategy) ---*/ @@ -434,7 +451,7 @@ 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); From 08a2f028eb9c6c75eb1a7eff29eec0b5c528668e Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 11:08:23 +0200 Subject: [PATCH 12/23] fix cfl for full multigrid --- SU2_CFD/include/solvers/CSolver.hpp | 6 +++ .../src/integration/CMultiGridIntegration.cpp | 54 ++++++++++++++----- 2 files changed, 47 insertions(+), 13 deletions(-) diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 942e2e25877..9d86e6afc5f 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -389,6 +389,12 @@ class CSolver { */ inline su2double GetAvg_CFL_Local(void) const { return Avg_CFL_Local; } + /*! + * \brief Set the value of the average local CFL number. + * \param[in] val_cfl - Average CFL number. + */ + inline void SetAvg_CFL_Local(su2double 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 b873605108e..6a5471024f1 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; @@ -67,6 +68,24 @@ static su2double GetMGLevelCorrectionScale(unsigned short iMesh) { } } +static passivedouble GetWarmupCFLScale(const std::vector& cflScaling, + unsigned long innerIter, + unsigned long startupIter) { + if (startupIter == 0 || cflScaling.empty()) return 1.0; + + const passivedouble firstScale = max(passivedouble{1e-6}, SU2_TYPE::GetValue(cflScaling.front())); + if (innerIter < startupIter) { + if (cflScaling.size() > 1) { + const passivedouble secondScale = max(passivedouble{1e-6}, SU2_TYPE::GetValue(cflScaling[1])); + return firstScale * secondScale; + } + return firstScale; + } + + if (innerIter < 2 * startupIter) return firstScale; + return 1.0; +} + inline passivedouble ComputeLinSysResRMS(const CSolver* solver) { passivedouble result = 0; for (unsigned short iVar = 0; iVar < solver->GetnVar(); ++iVar) { @@ -186,20 +205,29 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, SU2_OMP_SAFE_GLOBAL_ACCESS(config[iZone]->SubtractFinestMesh();) - auto nMG = config[iZone]->GetnMGLevels(); - auto cfl_base = SU2_TYPE::GetValue(config[iZone]->GetCFL(FinestMesh)); - auto cflScaling = config[iZone]->GetMGOptions().MG_CflScaling; - cout << "level =" << nMG - FinestMesh << endl; - cout <<"finestmesh = " << FinestMesh << endl; - cout << "CFL = " << cfl_base << ", scaling factor = " << cflScaling[FinestMesh-1] << endl; - // now scale the cfl using the scaling factors by multiplying from the first factor up to the factor for the current level - for (unsigned short iMesh = 0; iMesh < FinestMesh-1; ++iMesh) { - cout << "scale = " << cflScaling[iMesh] << endl; - cfl_base *= cflScaling[iMesh]; + /*--- Full-MG warmup: seed the history CFL with the base CFL scaled by the + * current warmup phase, so Avg_CFL_Local matches what the solver uses. ---*/ + static std::map base_cfl_by_zone; + if (base_cfl_by_zone.find(iZone) == base_cfl_by_zone.end()) { + base_cfl_by_zone[iZone] = SU2_TYPE::GetValue(config[iZone]->GetCFL(MESH_0)); } - cout << "new CFL = " << cfl_base << endl; - // set the new base cfl for MESH_0 to the scaled value: - SU2_OMP_SAFE_GLOBAL_ACCESS(config[iZone]->SetCFL(MESH_0, cfl_base);) + + const auto& cflScaling = config[iZone]->GetMGOptions().MG_CflScaling; + const unsigned long innerIter = config[iZone]->GetInnerIter(); + const passivedouble cfl_base = base_cfl_by_zone[iZone]; + const passivedouble warmupScale = GetWarmupCFLScale(cflScaling, innerIter, startup_iter); + const passivedouble warmupCFL = cfl_base * warmupScale; + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { + config[iZone]->SetCFL(MESH_0, warmupCFL); + CGeometry* geo_c = geometry[iZone][iInst][MESH_0]; + CSolver* sol_c = solver_container[iZone][iInst][MESH_0][Solver_Position]; + sol_c->SetAvg_CFL_Local(warmupCFL); + for (auto iPoint = 0ul; iPoint < geo_c->GetnPoint(); iPoint++) + sol_c->GetNodes()->SetLocalCFL(iPoint, warmupCFL); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS } From aac68022e0f83abd2759e03d96304e22394d7e30 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 11:33:29 +0200 Subject: [PATCH 13/23] fix cfl for full multigrid --- .../src/integration/CMultiGridIntegration.cpp | 70 ++++++------------- 1 file changed, 22 insertions(+), 48 deletions(-) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 6a5471024f1..1697356a817 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -204,64 +204,33 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, config[iZone]); SU2_OMP_SAFE_GLOBAL_ACCESS(config[iZone]->SubtractFinestMesh();) - - /*--- Full-MG warmup: seed the history CFL with the base CFL scaled by the - * current warmup phase, so Avg_CFL_Local matches what the solver uses. ---*/ - static std::map base_cfl_by_zone; - if (base_cfl_by_zone.find(iZone) == base_cfl_by_zone.end()) { - base_cfl_by_zone[iZone] = SU2_TYPE::GetValue(config[iZone]->GetCFL(MESH_0)); - } - - const auto& cflScaling = config[iZone]->GetMGOptions().MG_CflScaling; - const unsigned long innerIter = config[iZone]->GetInnerIter(); - const passivedouble cfl_base = base_cfl_by_zone[iZone]; - const passivedouble warmupScale = GetWarmupCFLScale(cflScaling, innerIter, startup_iter); - const passivedouble warmupCFL = cfl_base * warmupScale; - - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS - { - config[iZone]->SetCFL(MESH_0, warmupCFL); - CGeometry* geo_c = geometry[iZone][iInst][MESH_0]; - CSolver* sol_c = solver_container[iZone][iInst][MESH_0][Solver_Position]; - sol_c->SetAvg_CFL_Local(warmupCFL); - for (auto iPoint = 0ul; iPoint < geo_c->GetnPoint(); iPoint++) - sol_c->GetNodes()->SetLocalCFL(iPoint, warmupCFL); - } - END_SU2_OMP_SAFE_GLOBAL_ACCESS - - } /*--- Set the current finest grid (full multigrid strategy) ---*/ FinestMesh = config[iZone]->GetFinestMesh(); - /*--- Perform the Full Approximation Scheme multigrid ---*/ - - MultiGrid_Cycle(geometry, solver_container, numerics_container, config, - FinestMesh, RecursiveParam, RunTime_EqSystem, iZone, iInst); - - /*--- 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); @@ -269,17 +238,22 @@ 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_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); + /*--- 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], From 56134c50b822aa6580b6ba7c76c82c5ffda154b3 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 12:06:14 +0200 Subject: [PATCH 14/23] fix cfl for full multigrid --- SU2_CFD/include/solvers/CSolver.hpp | 10 ++++++++++ SU2_CFD/src/integration/CMultiGridIntegration.cpp | 1 + 2 files changed, 11 insertions(+) diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 9d86e6afc5f..67b7a7be69b 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -395,6 +395,16 @@ class CSolver { */ inline void SetAvg_CFL_Local(su2double val_cfl) { Avg_CFL_Local = val_cfl; } + /*! + * \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 1697356a817..8a721bf58f9 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -243,6 +243,7 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, 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); From 335e2c2bfc9174d8bf6e8f107a3fa51d3aa8d486 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 14:26:24 +0200 Subject: [PATCH 15/23] extra scaling for warmup phaser --- SU2_CFD/src/integration/CMultiGridIntegration.cpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 8a721bf58f9..e8d5198d522 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -223,6 +223,21 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, if (cfl_base < EPS) cfl_base = SU2_TYPE::GetValue(config[iZone]->GetCFL(MESH_0)); + /*--- FMG warmup: all warmup phases (FinestMesh > 0) solve a full problem from + * initial conditions on their active level; coarse-correction levels also + * use non-restricted ICs at first activation. Reduce the base CFL by + * 0.5^nMGLevels uniformly across all levels during warmup. This exactly + * replicates running with a reduced CFL_NUMBER during warmup and then + * recovering the full CFL when FinestMesh = 0 (V-cycle phase). + * Example (nMGLevels=2, CFL_NUMBER=100, MG_CFL_SCALING=0.3): + * Warmup phases 1 & 2: base=25, CFL[1]=7.5, CFL[2]=2.25 + * V-cycle phase: base=100, CFL[1]=30, CFL[2]=9 ---*/ + if (FinestMesh > 0) { + passivedouble warmup_factor = 1.0; + for (unsigned short k = 0; k < nMGLevels; ++k) warmup_factor *= 0.01; + cfl_base *= warmup_factor; + } + const auto& cflScaling = config[iZone]->GetMGOptions().MG_CflScaling; passivedouble CFL_local = cfl_base; From 006ecd57b8f706925932686dd500dad9bb83f9cd Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 19:48:53 +0200 Subject: [PATCH 16/23] cleanup implicit lines --- Common/src/geometry/CMultiGridGeometry.cpp | 99 ------------------- .../src/integration/CMultiGridIntegration.cpp | 15 --- .../integration/CSingleGridIntegration.cpp | 17 ++++ 3 files changed, 17 insertions(+), 114 deletions(-) diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index feeb1e2c60a..4936ae4bdbb 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1485,34 +1485,6 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, if (lines.empty()) return; - 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: * ANISOTROPIC (default): Pair nodes at the SAME distance from wall on DIFFERENT lines. * Each coarse CV has 2 fine children (from adjacent lines). @@ -1607,16 +1579,6 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& 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); @@ -1677,29 +1639,6 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, /*--- 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); @@ -1707,25 +1646,6 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& 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); @@ -1784,25 +1704,6 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, 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; } } diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index e8d5198d522..8a721bf58f9 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -223,21 +223,6 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, if (cfl_base < EPS) cfl_base = SU2_TYPE::GetValue(config[iZone]->GetCFL(MESH_0)); - /*--- FMG warmup: all warmup phases (FinestMesh > 0) solve a full problem from - * initial conditions on their active level; coarse-correction levels also - * use non-restricted ICs at first activation. Reduce the base CFL by - * 0.5^nMGLevels uniformly across all levels during warmup. This exactly - * replicates running with a reduced CFL_NUMBER during warmup and then - * recovering the full CFL when FinestMesh = 0 (V-cycle phase). - * Example (nMGLevels=2, CFL_NUMBER=100, MG_CFL_SCALING=0.3): - * Warmup phases 1 & 2: base=25, CFL[1]=7.5, CFL[2]=2.25 - * V-cycle phase: base=100, CFL[1]=30, CFL[2]=9 ---*/ - if (FinestMesh > 0) { - passivedouble warmup_factor = 1.0; - for (unsigned short k = 0; k < nMGLevels; ++k) warmup_factor *= 0.01; - cfl_base *= warmup_factor; - } - const auto& cflScaling = config[iZone]->GetMGOptions().MG_CflScaling; passivedouble CFL_local = cfl_base; diff --git a/SU2_CFD/src/integration/CSingleGridIntegration.cpp b/SU2_CFD/src/integration/CSingleGridIntegration.cpp index fab97842048..2806873d186 100644 --- a/SU2_CFD/src/integration/CSingleGridIntegration.cpp +++ b/SU2_CFD/src/integration/CSingleGridIntegration.cpp @@ -49,6 +49,23 @@ void CSingleGridIntegration::SingleGrid_Iteration(CGeometry ****geometry, CSolve CGeometry* geometry_fine = geometry[iZone][iInst][FinestMesh]; CSolver** solvers_fine = solver_container[iZone][iInst][FinestMesh]; + if (RunTime_EqSystem == RUNTIME_TURB_SYS) { + /*--- config[iZone]->GetCFL(FinestMesh) is already scaled for the active MG level + * by CMultiGridIntegration's per-level CFL rebuild (chain of MG_CFL_SCALING + * factors down to FinestMesh), including during FMG warmup. Do not re-apply + * that scaling here or it gets squared during warmup. ---*/ + const su2double turbReduction = SU2_TYPE::GetValue(config[iZone]->GetCFLRedCoeff_Turb()); + const su2double turbCFL = SU2_TYPE::GetValue(config[iZone]->GetCFL(FinestMesh)) * turbReduction; + auto* turbSolver = solvers_fine[Solver_Position]; + + SU2_OMP_SAFE_GLOBAL_ACCESS(turbSolver->SetCFL_Local_Stats(turbCFL);) + SU2_OMP_FOR_STAT(roundUpDiv(geometry_fine->GetnPoint(), omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < geometry_fine->GetnPoint(); ++iPoint) { + turbSolver->GetNodes()->SetLocalCFL(iPoint, turbCFL); + } + END_SU2_OMP_FOR + } + /*--- Preprocessing ---*/ solvers_fine[Solver_Position]->Preprocessing(geometry_fine, solvers_fine, config[iZone], From 1b84b72d76b90e68eeb19e6584fddb91dd07e1e2 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 19:59:05 +0200 Subject: [PATCH 17/23] cleanup implicit lines with claude --- Common/src/geometry/CMultiGridGeometry.cpp | 268 ++++++++------------- 1 file changed, 94 insertions(+), 174 deletions(-) diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 4936ae4bdbb..9749e1c76d6 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1393,7 +1393,6 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, 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 @@ -1432,7 +1431,10 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, norm_prev = sqrt(norm_prev); for (unsigned short d = 0; d < nDim; ++d) prev_dir[d] /= norm_prev; - /*--- Build the implicit line by following the best-aligned interior neighbor ---*/ + /*--- Build the implicit line by following the best-aligned interior neighbor. + * Nodes already on this line are excluded from candidacy so the walk + * cannot fold back on itself (the angle threshold alone would only make + * that unlikely, not impossible, on distorted meshes). ---*/ vector L; L.push_back(iPoint); auto current = iPoint; @@ -1442,10 +1444,10 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, unsigned long best_neighbor = ULONG_MAX; for (auto jPoint : fine_grid->nodes->GetPoints(current)) { - if (jPoint == current) continue; if (!fine_grid->nodes->GetDomain(jPoint)) continue; if (fine_grid->nodes->GetBoundary(jPoint)) continue; if (fine_grid->nodes->GetAgglomerate(jPoint)) continue; + if (find(L.begin(), L.end(), jPoint) != L.end()) continue; /*--- Compute normalized direction to candidate ---*/ su2double vec[MAXNDIM] = {0.0}; @@ -1485,198 +1487,115 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, if (lines.empty()) return; - /*--- 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. + /*--- Agglomeration strategy: at each "position" (distance from the wall along a + * line) every line contributes a block of nBlock consecutive nodes; two lines + * are paired through a mesh-adjacency search of their block anchors and the + * combined blocks become the children of one coarse CV. * - * 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. - ---*/ + * ANISOTROPIC (default, nBlock=1): pair single nodes at the SAME distance from + * the wall on two DIFFERENT lines. Each coarse CV has 2 fine children. + * Reduces the mesh by a factor ~2 normal to the wall, preserves resolution + * along the wall. + * + * ISOTROPIC (nBlock=2): pair 2-node blocks (2 positions x 2 lines) into one + * coarse CV with 4 fine children. Reduces the mesh uniformly by a factor + * ~4 in all directions. ---*/ + const unsigned long nBlock = ISOTROPIC ? 2 : 1; vector reserved(nPointFine, 0); - unsigned position_idx = 0; + unsigned long position_idx = 0; - while (true) { - bool any_work = false; - vector line_processed(lines.size(), 0); + /*--- Fine-grid nodes line `li` contributes at the current position, or empty if + the line is too short to reach that far. ---*/ + auto LineBlock = [&](unsigned long li) -> vector { + const auto& L = lines[li]; + const unsigned long first = 1 + nBlock * position_idx; + if (L.size() < first + nBlock) return {}; + return vector(L.begin() + first, L.begin() + first + nBlock); + }; - /*--- Build list of active lines (have nodes at current position) ---*/ + while (true) { + /*--- Cache each line's block for this position and collect the active ones. ---*/ + vector> block(lines.size()); 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; - 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 - } - active_lines.push_back(li); + block[li] = LineBlock(li); + if (!block[li].empty()) active_lines.push_back(li); } + if (active_lines.empty()) break; - if (ISOTROPIC) { - /*--- 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; - - 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; + vector line_processed(lines.size(), 0); + bool any_work = false; - /*--- 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; - } + for (auto li1 : active_lines) { + if (line_processed[li1]) 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); + const auto anchor = block[li1].front(); + if (fine_grid->nodes->GetAgglomerate(anchor) || reserved[anchor]) continue; - Index_CoarseCV++; - line_processed[li1] = line_processed[li2_best] = 1; - any_work = true; - } - } else { - /*--- 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; - } + /*--- Find an unprocessed neighboring line: one whose block anchor is a + mesh-neighbor of ours at the same position. ---*/ + unsigned long li2 = ULONG_MAX; + for (auto neighbor_point : fine_grid->nodes->GetPoints(anchor)) { + for (auto candidate : active_lines) { + if (candidate == li1 || line_processed[candidate]) continue; + if (block[candidate].front() == neighbor_point) { + li2 = candidate; + break; } - if (li2_best != std::numeric_limits::max()) break; } + if (li2 != ULONG_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; + if (li2 == ULONG_MAX) { + if (DEBUG_OUTPUT && position_idx < 3) { + cout << " Line " << li1 << " at position " << position_idx << " (node " << anchor + << ") has NO neighbor line!" << endl; } - - 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; - - /*--- 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_best] = 1; - any_work = true; + continue; } - } - position_idx++; - if (!any_work) break; + /*--- Assemble and validate the coarse CV's children together. ---*/ + auto group = block[li1]; + group.insert(group.end(), block[li2].begin(), block[li2].end()); - /*--- Check if any line still has positions available ---*/ - bool any_more = false; - if (ISOTROPIC) { - for (const auto& L : lines) { - if (L.size() > 1 + 2 * position_idx + 1) { - any_more = true; + bool valid = true; + for (auto p : group) { + if (fine_grid->nodes->GetAgglomerate(p) || reserved[p] || !GeometricalCheck(p, fine_grid, config)) { + valid = false; break; } } - } else { - for (const auto& L : lines) { - if (L.size() > 1 + position_idx) { - any_more = true; - break; - } + if (!valid) continue; // one of the nodes wasn't ready; li2 may still pair elsewhere. + + /*--- Guard against the same fine point appearing in both blocks + (can happen if two lines' walks overlap in space). ---*/ + bool duplicate = false; + for (size_t i = 0; !duplicate && i + 1 < group.size(); ++i) + for (size_t j = i + 1; j < group.size(); ++j) + if (group[i] == group[j]) duplicate = true; + + if (duplicate) { + line_processed[li1] = line_processed[li2] = 1; + continue; } + + /*--- Create the coarse CV from the combined blocks. ---*/ + for (unsigned short c = 0; c < group.size(); ++c) { + fine_grid->nodes->SetParent_CV(group[c], Index_CoarseCV); + nodes->SetChildren_CV(Index_CoarseCV, c, group[c]); + reserved[group[c]] = 1; + MGQueue_InnerCV.RemoveCV(group[c]); + } + nodes->SetnChildren_CV(Index_CoarseCV, static_cast(group.size())); + Index_CoarseCV++; + + line_processed[li1] = line_processed[li2] = 1; + any_work = true; } - if (!any_more) break; + + if (!any_work) break; + position_idx++; } /*--- Count how many CVs and nodes were created ---*/ @@ -1696,7 +1615,7 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, if (reserved[i]) nNodes_claimed++; } - if (rank == MASTER_NODE) { + if (DEBUG_OUTPUT) { 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; @@ -1707,14 +1626,15 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, } } - /*--- Verify all claimed nodes are properly marked as agglomerated ---*/ + /*--- Verify all claimed nodes are properly marked as agglomerated (SetParent_CV should + guarantee this; a mismatch would indicate a bookkeeping bug above). ---*/ 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) { + if (mismatches > 0 && DEBUG_OUTPUT) { cout << " WARNING: " << mismatches << " nodes marked as reserved but not agglomerated!" << endl; } } From 2c1edf08cc96b8039eebedd2903a7ec184cd51aa Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 20:00:29 +0200 Subject: [PATCH 18/23] cleanup implicit lines with claude --- .../src/integration/CMultiGridIntegration.cpp | 27 ------------------- 1 file changed, 27 deletions(-) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 8a721bf58f9..57c2b32cc7a 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -59,33 +59,6 @@ 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; - } -} - -static passivedouble GetWarmupCFLScale(const std::vector& cflScaling, - unsigned long innerIter, - unsigned long startupIter) { - if (startupIter == 0 || cflScaling.empty()) return 1.0; - - const passivedouble firstScale = max(passivedouble{1e-6}, SU2_TYPE::GetValue(cflScaling.front())); - if (innerIter < startupIter) { - if (cflScaling.size() > 1) { - const passivedouble secondScale = max(passivedouble{1e-6}, SU2_TYPE::GetValue(cflScaling[1])); - return firstScale * secondScale; - } - return firstScale; - } - - if (innerIter < 2 * startupIter) return firstScale; - return 1.0; -} - inline passivedouble ComputeLinSysResRMS(const CSolver* solver) { passivedouble result = 0; for (unsigned short iVar = 0; iVar < solver->GetnVar(); ++iVar) { From 38b9c11b75a1d7c79820b3adfe091bf00c43c23b Mon Sep 17 00:00:00 2001 From: Nijso Date: Wed, 5 Aug 2026 20:34:48 +0200 Subject: [PATCH 19/23] Apply suggestion from @bigfooted --- SU2_CFD/src/integration/CMultiGridIntegration.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 57c2b32cc7a..e34fd48149d 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -827,7 +827,7 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet 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]) -Residual_Fine[iVar] = 0.0; + Residual_Fine[iVar] = 0.0; su2double correction = factor * Residual_Fine[iVar]; Solution_Fine[iVar] += correction; From e9f8b7cd71cedd9157160edec5b6690498f758ea Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 20:44:01 +0200 Subject: [PATCH 20/23] remove debugging output --- Common/src/geometry/CMultiGridGeometry.cpp | 106 +----------------- .../src/integration/CMultiGridIntegration.cpp | 4 +- .../integration/CSingleGridIntegration.cpp | 5 +- SU2_CFD/src/iteration/CFluidIteration.cpp | 2 +- 4 files changed, 6 insertions(+), 111 deletions(-) diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 9749e1c76d6..7cea6d7d1c1 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -431,46 +431,6 @@ 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. ---*/ @@ -544,67 +504,6 @@ 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(); @@ -1431,10 +1330,7 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, norm_prev = sqrt(norm_prev); for (unsigned short d = 0; d < nDim; ++d) prev_dir[d] /= norm_prev; - /*--- Build the implicit line by following the best-aligned interior neighbor. - * Nodes already on this line are excluded from candidacy so the walk - * cannot fold back on itself (the angle threshold alone would only make - * that unlikely, not impossible, on distorted meshes). ---*/ + /*--- Build the implicit line by following the best-aligned interior neighbor. ---*/ vector L; L.push_back(iPoint); auto current = iPoint; diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index 57c2b32cc7a..f49674927cf 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -425,7 +425,6 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, solver_coarse->Preprocessing(geometry_coarse, solver_container_coarse, config, iMesh+1, NO_RK_ITER, RunTime_EqSystem, false); - 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) ---*/ @@ -452,6 +451,9 @@ 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); const auto& mgOpts = config->GetMGOptions(); diff --git a/SU2_CFD/src/integration/CSingleGridIntegration.cpp b/SU2_CFD/src/integration/CSingleGridIntegration.cpp index 2806873d186..0f051fdc2b6 100644 --- a/SU2_CFD/src/integration/CSingleGridIntegration.cpp +++ b/SU2_CFD/src/integration/CSingleGridIntegration.cpp @@ -50,10 +50,7 @@ void CSingleGridIntegration::SingleGrid_Iteration(CGeometry ****geometry, CSolve CSolver** solvers_fine = solver_container[iZone][iInst][FinestMesh]; if (RunTime_EqSystem == RUNTIME_TURB_SYS) { - /*--- config[iZone]->GetCFL(FinestMesh) is already scaled for the active MG level - * by CMultiGridIntegration's per-level CFL rebuild (chain of MG_CFL_SCALING - * factors down to FinestMesh), including during FMG warmup. Do not re-apply - * that scaling here or it gets squared during warmup. ---*/ + /*--- CFL scaling of turbulence during the warmup phase if FMG. ---*/ const su2double turbReduction = SU2_TYPE::GetValue(config[iZone]->GetCFLRedCoeff_Turb()); const su2double turbCFL = SU2_TYPE::GetValue(config[iZone]->GetCFL(FinestMesh)) * turbReduction; auto* turbSolver = solvers_fine[Solver_Position]; diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index 264effbdf7a..a79a9f1f8f4 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -98,7 +98,7 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe 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); + RUNTIME_TURB_SYS, val_iZone, val_iInst); } if (config[val_iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) { From 05e68d3dcde4c638fbde04ac3b33e5e36449b616 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 5 Aug 2026 22:22:50 +0200 Subject: [PATCH 21/23] remove unused variable --- Common/src/geometry/CMultiGridGeometry.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 7cea6d7d1c1..871d6fa8197 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -312,11 +312,8 @@ 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. ---*/ From f7677dcca4a71b16c61c4a8f91e108f68beb2933 Mon Sep 17 00:00:00 2001 From: Nijso Date: Wed, 5 Aug 2026 23:04:22 +0200 Subject: [PATCH 22/23] Potential fix for pull request finding 'CodeQL / Comparison of narrow type with wide type in loop condition' Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com> --- Common/src/geometry/CMultiGridGeometry.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 871d6fa8197..b97d919385b 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1474,7 +1474,7 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, } /*--- Create the coarse CV from the combined blocks. ---*/ - for (unsigned short c = 0; c < group.size(); ++c) { + for (size_t c = 0; c < group.size(); ++c) { fine_grid->nodes->SetParent_CV(group[c], Index_CoarseCV); nodes->SetChildren_CV(Index_CoarseCV, c, group[c]); reserved[group[c]] = 1; From ab003cb15ecf1305c40f274d2d233553ba4b1561 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Thu, 6 Aug 2026 13:47:09 +0200 Subject: [PATCH 23/23] update euler tests --- TestCases/euler/CRM/inv_CRM_JST.cfg | 16 ++++++++++------ TestCases/euler/channel/inv_channel.cfg | 9 ++++++--- TestCases/euler/naca0012/inv_NACA0012.cfg | 16 ++++++++++------ TestCases/euler/oneram6/inv_ONERAM6.cfg | 8 ++++---- 4 files changed, 30 insertions(+), 19 deletions(-) diff --git a/TestCases/euler/CRM/inv_CRM_JST.cfg b/TestCases/euler/CRM/inv_CRM_JST.cfg index 63c749021f2..60080a0254b 100644 --- a/TestCases/euler/CRM/inv_CRM_JST.cfg +++ b/TestCases/euler/CRM/inv_CRM_JST.cfg @@ -43,11 +43,11 @@ MARKER_MONITORING= ( fuselage , Wing , HTP ) % ------------- COMMON PARAMETERS TO DEFINE THE NUMERICAL METHOD --------------% % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -CFL_NUMBER= 5.0 +CFL_NUMBER= 100.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) -EXT_ITER= 99999 +ITER= 1000 LINEAR_SOLVER= FGMRES LINEAR_SOLVER_ERROR= 1E-1 LINEAR_SOLVER_ITER= 3 @@ -55,12 +55,16 @@ LINEAR_SOLVER_ITER= 3 % -------------------------- MULTIGRID PARAMETERS -----------------------------% % MGLEVEL= 3 +MG_MIN_MESHSIZE= 100 MGCYCLE= V_CYCLE -MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) -MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_SMOOTH_OUTPUT= YES +MG_SMOOTH_EARLY_EXIT= YES +MG_PRE_SMOOTH= ( 5, 5, 5, 5 ) +MG_POST_SMOOTH= ( 5, 5, 5, 5 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) -MG_DAMP_RESTRICTION= 0.5 -MG_DAMP_PROLONGATION= 0.5 +MG_DAMP_RESTRICTION= 0.75 +MG_DAMP_PROLONGATION= 0.75 +MG_CFL_SCALING= 0.5, 0.5, 0.5 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/euler/channel/inv_channel.cfg b/TestCases/euler/channel/inv_channel.cfg index 5161a04413d..8145521b147 100644 --- a/TestCases/euler/channel/inv_channel.cfg +++ b/TestCases/euler/channel/inv_channel.cfg @@ -47,11 +47,11 @@ MARKER_MONITORING= ( upper_wall, lower_wall ) % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % NUM_METHOD_GRAD= GREEN_GAUSS -CFL_NUMBER= 6.0 +CFL_NUMBER= 100.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) -EXT_ITER= 999999 +ITER= 1000 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % @@ -63,12 +63,15 @@ LINEAR_SOLVER_ITER= 3 % -------------------------- MULTIGRID PARAMETERS -----------------------------% % MGLEVEL= 3 -MGCYCLE= V_CYCLE +MGCYCLE= W_CYCLE +MG_SMOOTH_OUTPUT= NO +MG_SMOOTH_EARLY_EXIT= YES MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) MG_POST_SMOOTH= ( 4, 4, 4, 4 ) MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) MG_DAMP_RESTRICTION= 0.5 MG_DAMP_PROLONGATION= 0.5 +MG_CFL_SCALING= 0.5, 0.5, 0.5 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/euler/naca0012/inv_NACA0012.cfg b/TestCases/euler/naca0012/inv_NACA0012.cfg index ca235bb07a3..814649dbc30 100644 --- a/TestCases/euler/naca0012/inv_NACA0012.cfg +++ b/TestCases/euler/naca0012/inv_NACA0012.cfg @@ -46,7 +46,7 @@ MARKER_DESIGNING = ( airfoil ) % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES OBJECTIVE_FUNCTION= DRAG -CFL_NUMBER= 5.0 +CFL_NUMBER= 100.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) ITER= 250 @@ -62,11 +62,15 @@ LINEAR_SOLVER_ITER= 5 % MGLEVEL= 3 MGCYCLE= W_CYCLE -MG_PRE_SMOOTH= ( 1, 2, 3, 3 ) -MG_POST_SMOOTH= ( 0, 0, 0, 0 ) -MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 ) -MG_DAMP_RESTRICTION= 1.0 -MG_DAMP_PROLONGATION= 1.0 +MG_SMOOTH_OUTPUT= YES +MG_SMOOTH_EARLY_EXIT= YES + +MG_PRE_SMOOTH= ( 5, 5, 5, 5 ) +MG_POST_SMOOTH= ( 5, 5, 5, 5 ) +MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) +MG_DAMP_RESTRICTION= 0.5 +MG_DAMP_PROLONGATION= 0.5 +MG_CFL_SCALING= 0.5, 0.5, 0.5 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/euler/oneram6/inv_ONERAM6.cfg b/TestCases/euler/oneram6/inv_ONERAM6.cfg index 6517b957871..59c19e14a53 100644 --- a/TestCases/euler/oneram6/inv_ONERAM6.cfg +++ b/TestCases/euler/oneram6/inv_ONERAM6.cfg @@ -49,15 +49,15 @@ MARKER_DESIGNING = ( WING ) % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES OBJECTIVE_FUNCTION= DRAG -CFL_NUMBER= 5.0 +CFL_NUMBER= 100.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) ITER= 99999 LINEAR_SOLVER= FGMRES LINEAR_SOLVER_PREC= LU_SGS -LINEAR_SOLVER_ERROR= 1E-6 -LINEAR_SOLVER_ITER= 2 +LINEAR_SOLVER_ERROR= 1E-1 +LINEAR_SOLVER_ITER= 5 % ----------------------- SLOPE LIMITER DEFINITION ----------------------------% % @@ -92,7 +92,7 @@ TIME_DISCRE_ADJFLOW= EULER_IMPLICIT % --------------------------- CONVERGENCE PARAMETERS --------------------------& % CONV_RESIDUAL_MINVAL= -12 -CONV_STARTITER= 25 +CONV_STARTITER= 10 CONV_CAUCHY_ELEMS= 100 CONV_CAUCHY_EPS= 1E-10