diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 2dcd241f9de..5bda4780d23 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1126,6 +1126,8 @@ 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. */ + 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 a9ca92ab145..2d68a229b27 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2069,6 +2069,10 @@ 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_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 684da742b13..feeb1e2c60a 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,46 @@ 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 +544,67 @@ 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 +763,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); @@ -1286,8 +1388,12 @@ 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(); + 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 @@ -1381,55 +1487,103 @@ 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; + } } - /*--- 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; + 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; - const auto idx2 = 1 + 2 * pair_idx + 1; - if (L.size() <= idx2) continue; // no pair at this stage - const auto pW = fine_grid->nodes->GetParent_CV(L[0]); - parent_to_lines[pW].push_back(li); + 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); } - 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; + 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& L2 = lines[li2]; - const auto idx1 = 1 + 2 * pair_idx; + const auto idx1 = 1 + 2 * position_idx; const auto idx2 = idx1 + 1; - if (L1.size() <= idx2 || L2.size() <= idx2) continue; + 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(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; + 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) || @@ -1438,11 +1592,11 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, /*--- 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; + line_processed[li1] = line_processed[li2_best] = 1; continue; } - /*--- Create 4-child coarse CV ---*/ + /*--- 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); @@ -1453,6 +1607,16 @@ 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); @@ -1460,52 +1624,196 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV, 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; + 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; - /*--- 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); + 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; + } - Index_CoarseCV++; - any_work = true; + 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; + } } - 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; } + + /*--- 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..a14ad116f03 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 { @@ -53,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) { @@ -154,12 +168,14 @@ 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 ))) { + 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], @@ -197,10 +213,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); @@ -283,8 +299,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(); @@ -411,8 +432,9 @@ 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); + 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); @@ -437,8 +459,6 @@ 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(); @@ -814,7 +834,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; @@ -852,21 +872,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 diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index a79a9f1f8f4..264effbdf7a 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) { diff --git a/config_template.cfg b/config_template.cfg index bc6be98103a..45a822755c9 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1723,6 +1723,15 @@ 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 +% +% 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 -----------------------------% %