Multigrid with 3D implicit line agglomeration by bigfooted · Pull Request #2863 · su2code/SU2

bigfooted · 2026-08-05T18:04:36Z

Proposed Changes

fix implicit line agglomeration

When MG_IMPLICIT_LINE_AGGLOMERATION=YES, we try to create structured coarser meshes by first agglomerating the boundary, then use a paving method, also known as an advancing front method. We advance the entire front layer by layer, so all viscous walls advance one layer at a time, until the end, until an edge direction deviates more than 30 degrees or until a collision. When paving for a single cell fails, paving for the entire marker stops.

After agglomeration of the viscous walls, the other boundaries are agglomerated in a structured way if their first cell aspect ratio is higher than 2. For structured quadrilateral meshes with inflation layers this method leads to very high quality multigrids.
With MPI, the layers try to continue through the mpi interface. When the mpi interface created by parmetis is not very 'straight', this method often terminates and quality deteriorates. However, even the fallback to regular agglomeration still leads to high quality meshes.

Related Work

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

I am submitting my contribution to the develop branch.
My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
I have added a test case that demonstrates my contribution, if necessary.
I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

remove multigrid turbulence Co-authored-by: Nijso <nijso@hotmail.com>

Co-authored-by: Nijso <nijso@hotmail.com>

…feature_fullmg

… type with wide type in loop condition' Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>

bigfooted · 2026-09-06T19:19:40Z

screenshot of the first level agglomeration for the turbulent flat plate. the black lines is the original mesh, the red lines is the new agglomerated mesh. Note that boundary nodes are not agglomerated with interior nodes.

bigfooted · 2026-09-06T19:24:51Z

3D bump in channel, vertical slice through the domain:

bigfooted · 2026-09-06T19:30:57Z

This image is a close-up of the slice through the 3D bump in channel. The orange line shows the location of the MPI interface. Left and right of the bump the wall becomes a symmetry plane, and at this interface there is no horizontal agglomeration.

bigfooted · 2026-09-06T19:35:51Z

note that most of the agglomerations where no horizontal agglomeration occurs is due to some MPI collision. without MPI, the agglomeration looks better:

This shows that improving parmetis MPI agglomeration can result in better multigrid agglomeration

…_fullmg

Co-authored-by: Nijso <nijso@hotmail.com>

bigfooted · 2026-09-07T18:02:39Z

Code is still a bit bloated, I will do a second pass...
I will also modify some testcase setups so they have settings that converge well. Some testcases you can increase CFL from 10 -> 1000 without any problem.

+    groups.push_back({si});
+  }
+
+  const unsigned nRounds = (max_group <= 1) ? 0 : ((max_group <= 2) ? 1 : 2);


bigfooted · 2026-09-08T20:12:25Z

RAE2822 airfoil:

and a zoom:

bigfooted · 2026-09-08T20:33:44Z

RAE2822 case
convergence with implicit lines: 185 iterations, 19 s.

single grid: 740 iterations, 29 s. CFL had to be reduced from 1000 -> 100 to converge it. And it still has some weird convergence peaks at the end

+  for (unsigned long layer = 1;; ++layer) {
+    /*--- Every rank runs the same number of rounds, each ends in a collective handover
+     *    exchange. ---*/
+    int aliveLocal = 0;
+    for (unsigned long f = 0; f < fronts.size(); ++f) aliveLocal |= fronts[f].alive;
+    int aliveGlobal = 0;
+    SU2_MPI::Allreduce(&aliveLocal, &aliveGlobal, 1, MPI_INT, MPI_MAX, SU2_MPI::GetComm());
+    if (aliveGlobal == 0) break;
+
+    for (auto& F : fronts) F.failed = F.keepLocal = F.handTag = 0;
+    for (const auto& b : bids) bidIdx[b.node] = NOBID;
+    bids.clear();
+    bidOwner.clear();
+
+    /*--- (a) Every alive front proposes a successor for each of its nodes. A front that cannot fill
+     *    a whole layer proposes nothing and retires this round. ---*/
+    for (unsigned long f = 0; f < fronts.size(); ++f) {
+      if (!fronts[f].alive) continue;
+      fronts[f].prop.clear();
+      fronts[f].handTo.clear();
+
+      for (auto n : fronts[f].nodes) {
+        const auto c = bestSuccessor(n, fronts[f].dir.data());
+
+        /*--- Nothing free here but the stack continues across the interface. ---*/
+        if ((c.node == NO_POINT) && (c.halo != NO_POINT)) {
+          fronts[f].handTo.push_back(c.halo);
          continue;
+        }
+        /*--- No successor at all: a boundary, a partition, another front, or unusable mesh. ---*/
+        if (c.node == NO_POINT) {
+          markFail(f);
+          break;
+        }

-        /*--- 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;
+        CStep s{c.node, n, fine_grid->nodes->GetGlobalIndex(n), c.dot, c.len, {}};
+        for (unsigned short d = 0; d < nDim; ++d) s.dir[d] = c.dir[d];
+        fronts[f].prop.push_back(s);
+      }
+
+      /*--- An interface can cut a footprint. All of it crossing hands the stack over intact, part of
+       *    it crossing splits the footprint into a piece that stays and a piece handed across. ---*/
+      if (fronts[f].failed) {
+        fronts[f].prop.clear();
+        fronts[f].handTo.clear();
+      } else if (!fronts[f].handTo.empty()) {
+        /*--- fronts[f].prop is built in the order of fronts[f].nodes, so this is the piece that
+         *    stays, in the same order as the layer it proposes. ---*/
+        vector<unsigned long> narrow;
+        for (const auto& s : fronts[f].prop) narrow.push_back(s.from);
+
+        /*--- A cut can leave the local piece disconnected, which is not a layer. Drop it and hand
+         *    over the rest. ---*/
+        if (!narrow.empty() && !IsConnectedLayer(fine_grid, narrow)) {
+          narrow.clear();
+        }
+
+        fronts[f].handTag = TagOfSet(fine_grid, fronts[f].handTo);
+
+        if (narrow.empty()) {
+          fronts[f].prop.clear();
+        } else {
+          /*--- Close the coarse CV open on the wide footprint before narrowing, so that no CV holds
+           *    two layers of different shape. ---*/
+          fronts[f].nodes = narrow;
+          emit(f);
+          fronts[f].nBlock = BlockFor(maxAgglomSize, fronts[f].nodes);
+          fronts[f].tag = TagOfSet(fine_grid, fronts[f].nodes);
+          fronts[f].keepLocal = 1;
+        }
+      }
+    }
+
+    /*--- (b) Contention resolved from bids that were all collected before any was granted, so the
+     *    outcome does not depend on the order the fronts are visited in. ---*/
+    auto better = [](const CStep& a, const CStep& b) {
+      if (a.score != b.score) return a.score > b.score;
+      if (a.dist != b.dist) return a.dist < b.dist;
+      return a.key < b.key;
+    };
+
+    for (unsigned long f = 0; f < fronts.size(); ++f) {
+      if (!fronts[f].alive || fronts[f].failed) continue;
+      for (const auto& s : fronts[f].prop) {
+        /*--- A front that has lost a bid retires and places no more of its layer. ---*/
+        if (fronts[f].failed) break;
+
+        if (bidIdx[s.node] == NOBID) {
+          bidIdx[s.node] = static_cast<unsigned>(bids.size());
+          bids.push_back(s);
+          bidOwner.push_back(f);
          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);
+        const auto k = bidIdx[s.node];
+        const auto g = bidOwner[k];
+        /*--- Two nodes of the same front reaching for one successor is a pinch, the layer would come
+         *    out narrower than the front. ---*/
+        if (better(s, bids[k])) {
+          markFail(g);
+          bids[k] = s;
+          bidOwner[k] = f;
+        } else {
+          markFail(f);
+        }

-        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;
+        /*--- A head-on meeting stops both fronts, a glancing contact only costs the loser. ---*/
+        if ((g != f) && (GeometryToolbox::DotProduct(nDim, fronts[f].dir.data(), fronts[g].dir.data()) < 0.0)) {
+          markFail(f);
+          markFail(g);
+        }
      }
    }

-    /*--- 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);
+    /*--- (c) All-or-nothing acceptance: a front takes the whole layer or none of it and retires. A
+     *    bid only becomes a claim here. ---*/
+    for (unsigned long f = 0; f < fronts.size(); ++f) {
+      if (!fronts[f].alive) continue;
+
+      newLayer.clear();
+      if (!fronts[f].failed) {
+        /*--- Built in proposal order, so newLayer[k] is the successor of nodes[k]. ---*/
+        for (const auto& s : fronts[f].prop) {
+          const auto k = bidIdx[s.node];
+          if ((k != NOBID) && (bidOwner[k] == f)) newLayer.push_back(s.node);
+        }
+        /*--- Every bid of a front that was not marked failed must have been granted. ---*/
+        if ((newLayer.size() != fronts[f].nodes.size()) || !LayerIsIsomorphic(fine_grid, fronts[f].nodes, newLayer))
+          markFail(f);
+      }

-      Index_CoarseCV++;
-      any_work = true;
+      if (fronts[f].failed) {
+        fronts[f].alive = 0;
+        emit(f);
+        continue;
+      }
+
+      /*--- Turn the marching direction towards the mean of the steps just taken. ---*/
+      su2double mean[MAXNDIM] = {0.0};
+      for (const auto& s : fronts[f].prop)
+        for (unsigned short d = 0; d < nDim; ++d) mean[d] += s.dir[d];
+      const su2double meanNrm = GeometryToolbox::Norm(nDim, mean);
+      if (meanNrm > 0.0) {
+        su2double blended[MAXNDIM] = {0.0};
+        for (unsigned short d = 0; d < nDim; ++d)
+          blended[d] = (1.0 - DIR_BLEND) * fronts[f].dir[d] + DIR_BLEND * mean[d] / meanNrm;
+        const su2double bNrm = GeometryToolbox::Norm(nDim, blended);
+        if (bNrm > 0.0)
+          for (unsigned short d = 0; d < nDim; ++d) fronts[f].dir[d] = blended[d] / bNrm;
+      }
+
+      for (auto p : newLayer) {
+        claimed[p] = 1;
+      }
+      fronts[f].nodes = std::move(newLayer);
+      fronts[f].depth++;
+      ct[P_LAYERS]++;
+
+      fronts[f].pending.insert(fronts[f].pending.end(), fronts[f].nodes.begin(), fronts[f].nodes.end());
+      fronts[f].pendingLayers++;
+      if (fronts[f].pendingLayers >= fronts[f].nBlock) emit(f);
    }

-    pair_idx++;
-    if (!any_work) break;
+    /*--- (d) Hand stacks across partition interfaces. The footprint is sent to the owning rank,
+     *    packed against the receive marker and sent to the rank that marker receives from. ---*/
+    for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
+      if (!((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && (config->GetMarker_All_SendRecv(iMarker) > 0)))
+        continue;

-    /*--- 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;
+      const auto MarkerS = iMarker, MarkerR = iMarker + 1;
+      const auto send_to = config->GetMarker_All_SendRecv(MarkerS) - 1;
+      const auto receive_from = abs(config->GetMarker_All_SendRecv(MarkerR)) - 1;
+      const auto nVertexS = fine_grid->nVertex[MarkerS];
+      const auto nVertexR = fine_grid->nVertex[MarkerR];
+
+      /*--- Packed against the halo vertices. Tag and direction go separately, the AD MPI wrapper has
+       *    no byte type to send a struct. ---*/
+      vector<unsigned long> tagOut(nVertexR, 0), tagIn(nVertexS, 0);
+      vector<su2double> dirOut(nVertexR * nDim, 0.0), dirIn(nVertexS * nDim, 0.0);
+
+      for (auto& F : fronts)
+        for (auto p : F.handTo) {
+          if (haloMarker[p] != static_cast<int>(MarkerR)) continue;
+          const auto v = haloVertex[p];
+          /*--- Two fronts reaching for one node: the lower tag takes it. ---*/
+          if ((tagOut[v] != 0) && (tagOut[v] <= F.handTag)) continue;
+          tagOut[v] = F.handTag;
+          for (unsigned short d = 0; d < nDim; ++d) dirOut[v * nDim + d] = F.dir[d];
+        }
+
+      SU2_MPI::Sendrecv(tagOut.data(), nVertexR, MPI_UNSIGNED_LONG, receive_from, 2, tagIn.data(), nVertexS,
+                        MPI_UNSIGNED_LONG, send_to, 2, SU2_MPI::GetComm(), MPI_STATUS_IGNORE);
+      SU2_MPI::Sendrecv(dirOut.data(), nVertexR * nDim, MPI_DOUBLE, receive_from, 3, dirIn.data(), nVertexS * nDim,
+                        MPI_DOUBLE, send_to, 3, SU2_MPI::GetComm(), MPI_STATUS_IGNORE);
+
+      for (auto iVertex = 0ul; iVertex < nVertexS; iVertex++) {
+        if (tagIn[iVertex] == 0) continue;
+        inherited.push_back({tagIn[iVertex], fine_grid->vertex[MarkerS][iVertex]->GetNode(), {}});
+        for (unsigned short d = 0; d < nDim; ++d) inherited.back().dir[d] = dirIn[iVertex * nDim + d];
      }
    }
-    if (!any_more) break;
+
+    /*--- A front that handed its whole footprint over is finished here, one that handed over only a
+     *    piece keeps marching on what was left. ---*/
+    for (unsigned long f = 0; f < fronts.size(); ++f) {
+      if (fronts[f].handTo.empty()) continue;
+      fronts[f].handTo.clear();
+      if (fronts[f].keepLocal) continue;
+      fronts[f].alive = 0;
+      emit(f);
+    }
+
+    /*--- Adopt what the neighbours sent, tags ascending so arrival order cannot change the outcome.
+     *    A footprint whose nodes are not all free is dropped. ---*/
+    std::sort(inherited.begin(), inherited.end(),
+              [](const CInherited& a, const CInherited& b) { return a.tag < b.tag; });
+
+    for (size_t i = 0; i < inherited.size();) {
+      size_t j = i;
+      while ((j < inherited.size()) && (inherited[j].tag == inherited[i].tag)) ++j;
+
+      vector<unsigned long> layer0;
+      bool ok = true;
+      for (size_t k = i; k < j; ++k) {
+        const auto p = inherited[k].node;
+        if (claimed[p] || fine_grid->nodes->GetAgglomerate(p) || !GeometricalCheck(p, fine_grid, config)) ok = false;
+        layer0.push_back(p);
+      }
+      /*--- The footprint has to arrive whole and connected. ---*/
+      if (ok && !IsConnectedLayer(fine_grid, layer0)) ok = false;
+
+      if (ok) {
+        std::array<su2double, MAXNDIM> d0{};
+        for (unsigned short d = 0; d < nDim; ++d) d0[d] = inherited[i].dir[d];
+        /*--- An inherited layer is an interior one, so it opens an ordinary two-deep coarse CV. ---*/
+        const auto nf = addFront(layer0, d0, inherited[i].tag, BlockFor(maxAgglomSize, layer0));
+        for (auto p : layer0) {
+          claimed[p] = 1;
+        }
+        ct[P_LAYERS]++;
+        if (fronts[nf].pendingLayers >= fronts[nf].nBlock) emit(nf);
+      }
+      i = j;
+    }
+    inherited.clear();
+  }


bigfooted and others added 20 commits July 29, 2026 19:36

fix order

9a2d329

minor implicit line changes

06790fb

Merge branch 'develop' into fix_mg_implicit_lines

4cc6571

fix multigrid turbulence

3b9e9df

small flow solver update

7e1f670

small flow solver update

33f2eaf

Apply suggestions from code review

3e0d528

remove multigrid turbulence Co-authored-by: Nijso <nijso@hotmail.com>

Merge branch 'develop' into feature_fullmg

feccac0

Apply suggestions from code review

4922e7b

Co-authored-by: Nijso <nijso@hotmail.com>

precommit

5d80679

Apply suggestions from code review

7a0db10

Co-authored-by: Nijso <nijso@hotmail.com>

Apply suggestions from code review

aca34af

Co-authored-by: Nijso <nijso@hotmail.com>

fix cfl for full multigrid

5648a3c

fix cfl for full multigrid

08a2f02

fix cfl for full multigrid

aac6802

fix cfl for full multigrid

56134c5

extra scaling for warmup phaser

335e2c2

cleanup implicit lines

006ecd5

cleanup implicit lines with claude

1b84b72

cleanup implicit lines with claude

2c1edf0

bigfooted added changelog:fix PR: draft labels Aug 5, 2026

bigfooted commented Aug 5, 2026

View reviewed changes

Comment thread SU2_CFD/src/integration/CMultiGridIntegration.cpp Outdated

bigfooted added 3 commits August 5, 2026 20:34

Apply suggestion from @bigfooted

38b9c11

remove debugging output

e9f8b7c

Merge branch 'feature_fullmg' of https://github.com/su2code/su2 into …

98da9bc

…feature_fullmg

github-advanced-security Bot found potential problems Aug 5, 2026

View reviewed changes

Comment thread Common/src/geometry/CMultiGridGeometry.cpp Fixed

Comment thread Common/src/geometry/CMultiGridGeometry.cpp Fixed

bigfooted and others added 3 commits August 5, 2026 22:22

remove unused variable

05e68d3

Potential fix for pull request finding 'CodeQL / Comparison of narrow…

f7677dc

… type with wide type in loop condition' Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>

update euler tests

ab003cb

bigfooted added 3 commits September 6, 2026 19:00

cleanup diagnostics

625317b

cleanup config option

1251699

Merge branch 'develop' into feature_fullmg

1cf2dc4

bigfooted added 5 commits September 6, 2026 22:06

cleanup comments

6008bd3

Merge branch 'develop' of https://github.com/su2code/su2 into feature…

ccf853d

…_fullmg

properly write the implicit line agglomeration info

5de0a00

fix bug with cfl scaling with turbulence

3fc6736

Merge branch 'develop' of https://github.com/su2code/su2 into feature…

44bc799

…_fullmg

bigfooted removed the PR: draft label Sep 7, 2026

bigfooted changed the title ~~[WIP] Multigrid with 3D implicit line agglomeration~~ Multigrid with 3D implicit line agglomeration Sep 7, 2026