[WIP] Full multigrid with implicit line agglomeration

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. */
 };
 
 /*!

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);
 

Common/src/geometry/CMultiGridGeometry.cpp

@@ -659,11 +659,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 +1284,11 @@ 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 ---*/
 
   /*--- Collect implicit lines starting at viscous (no-slip) wall vertices only.
    *    Seeding from non-wall boundaries (farfield, inlet, outlet, symmetry) would
@@ -1326,7 +1327,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 ---*/
+      /*--- Build the implicit line by following the best-aligned interior neighbor. ---*/
       vector<unsigned long> L;
       L.push_back(iPoint);
       auto current = iPoint;
@@ -1336,10 +1337,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};
@@ -1379,133 +1380,154 @@ void CMultiGridGeometry::AgglomerateImplicitLines(unsigned long& Index_CoarseCV,
 
   if (lines.empty()) return;
 
-  if (rank == MASTER_NODE) {
-    cout << "Implicit line agglomeration: detected " << lines.size() << " lines." << 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: 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.
+   *
+   *    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<char> reserved(nPointFine, 0);
-  unsigned pair_idx = 0;
+  unsigned long position_idx = 0;
 
-  while (true) {
-    bool any_work = false;
+  /*--- 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<unsigned long> {
+    const auto& L = lines[li];
+    const unsigned long first = 1 + nBlock * position_idx;
+    if (L.size() < first + nBlock) return {};
+    return vector<unsigned long>(L.begin() + first, L.begin() + first + nBlock);
+  };
 
-    /*--- Build map: wall parent CV -> list of line indices ---*/
-    unordered_map<unsigned long, vector<unsigned long>> parent_to_lines;
-    parent_to_lines.reserve(lines.size());
+  while (true) {
+    /*--- Cache each line's block for this position and collect the active ones. ---*/
+    vector<vector<unsigned long>> block(lines.size());
+    vector<unsigned long> 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);
+      block[li] = LineBlock(li);
+      if (!block[li].empty()) active_lines.push_back(li);
     }
+    if (active_lines.empty()) break;
 
     vector<char> line_processed(lines.size(), 0);
+    bool any_work = false;
+
+    for (auto li1 : active_lines) {
+      if (line_processed[li1]) continue;
+
+      const auto anchor = block[li1].front();
+      if (fine_grid->nodes->GetAgglomerate(anchor) || reserved[anchor]) continue;
+
+      /*--- 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 != ULONG_MAX) break;
+      }
 
-    /*--- 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;
+      if (li2 == ULONG_MAX) {
+        if (DEBUG_OUTPUT && position_idx < 3) {
+          cout << "  Line " << li1 << " at position " << position_idx << " (node " << anchor
+               << ") has NO neighbor line!" << endl;
         }
+        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);
+      /*--- Assemble and validate the coarse CV's children together. ---*/
+      auto group = block[li1];
+      group.insert(group.end(), block[li2].begin(), block[li2].end());
 
-        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;
+      bool valid = true;
+      for (auto p : group) {
+        if (fine_grid->nodes->GetAgglomerate(p) || reserved[p] || !GeometricalCheck(p, fine_grid, config)) {
+          valid = false;
+          break;
+        }
       }
-    }
+      if (!valid) continue;  // one of the nodes wasn't ready; li2 may still pair elsewhere.
 
-    /*--- 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);
+      /*--- 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 (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;
+        MGQueue_InnerCV.RemoveCV(group[c]);
+      }
+      nodes->SetnChildren_CV(Index_CoarseCV, static_cast<unsigned short>(group.size()));
       Index_CoarseCV++;
+
+      line_processed[li1] = line_processed[li2] = 1;
       any_work = true;
     }
 
-    pair_idx++;
     if (!any_work) break;
+    position_idx++;
+  }
 
-    /*--- Check if any line still has pairs at the next stage ---*/
-    bool any_more = false;
-    for (const auto& L : lines) {
-      if (L.size() > 1 + 2 * pair_idx + 1) {
-        any_more = true;
-        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++;
     }
-    if (!any_more) break;
+  }
+
+  for (unsigned long i = 0; i < nPointFine; ++i) {
+    if (reserved[i]) nNodes_claimed++;
+  }
+
+  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;
+
+    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;
+    }
+  }
+
+  /*--- 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 && DEBUG_OUTPUT) {
+    cout << "  WARNING: " << mismatches << " nodes marked as reserved but not agglomerated!" << endl;
   }
 }

SU2_CFD/include/solvers/CSolver.hpp

@@ -389,6 +389,22 @@ 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 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.
    */