Multigrid with 3D implicit line agglomeration - #2863
Open
bigfooted wants to merge 69 commits into
Open
Conversation
remove multigrid turbulence Co-authored-by: Nijso <nijso@hotmail.com>
Co-authored-by: Nijso <nijso@hotmail.com>
Co-authored-by: Nijso <nijso@hotmail.com>
Co-authored-by: Nijso <nijso@hotmail.com>
bigfooted
commented
Aug 5, 2026
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
… type with wide type in loop condition' Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>
Contributor
Author
|
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. |
Contributor
Author
|
3D bump in channel, vertical slice through the domain: |
Contributor
Author
|
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. |
Contributor
Author
|
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 |
bigfooted
changed the title
bigfooted
commented
Sep 7, 2026
bigfooted
commented
Sep 7, 2026
bigfooted
commented
Sep 7, 2026
bigfooted
commented
Sep 7, 2026
bigfooted
commented
Sep 7, 2026
bigfooted
commented
Sep 7, 2026
bigfooted
commented
Sep 7, 2026
bigfooted
commented
Sep 7, 2026
Co-authored-by: Nijso <nijso@hotmail.com>
Contributor
Author
|
Code is still a bit bloated, I will do a second pass... |
| groups.push_back({si}); | ||
| } | ||
|
|
||
| const unsigned nRounds = (max_group <= 1) ? 0 : ((max_group <= 2) ? 1 : 2); |
Comment on lines
+2213
to
+2471
| for (unsigned long layer = 1;; ++layer) { | ||
| /*--- Every rank runs the same number of rounds: each ends in a collective handover exchange, so | ||
| * one dropping out early would hang the others. ---*/ | ||
| 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: it is retiring this round anyway, and letting its partial bids | ||
| * stand would let a dying front displace a healthy one out of nodes it can still use. ---*/ | ||
| 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; the split test below | ||
| * decides whether the whole layer goes over. ---*/ | ||
| 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; | ||
| } | ||
|
|
||
| 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; | ||
| 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); | ||
| } | ||
|
|
||
| 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; | ||
| /*--- An interface can cut a footprint. If all of it crosses, the stack is handed over intact | ||
| * and this front is finished. If only part crosses, the footprint is SPLIT: the piece whose | ||
| * successors are local marches on here, the rest is handed across, and both are renamed. ---*/ | ||
| 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, and phi still runs index for index between it and 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 in two disconnected halves - a square footprint cut | ||
| * diagonally does exactly that - and that is not a layer. Drop it and hand over the rest; | ||
| * the stack still survives on the far side instead of ending here. ---*/ | ||
| if (!narrow.empty() && !IsConnectedLayer(fine_grid, narrow)) { | ||
| narrow.clear(); | ||
| } | ||
|
|
||
| const auto a = L1[idx1], b = L1[idx2]; | ||
| const auto c = L2[idx1], d = L2[idx2]; | ||
| fronts[f].handTag = TagOfSet(fine_grid, fronts[f].handTo); | ||
|
|
||
| if (narrow.empty()) { | ||
| fronts[f].prop.clear(); | ||
| } else { | ||
| /*--- Close the coarse CV that is open on the WIDE footprint before narrowing, so that no CV | ||
| * ever ends up holding 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; | ||
| } | ||
| } | ||
| } | ||
|
|
||
| /*--- 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)) | ||
| /*--- (b) Contention resolved from bids that were all collected before any was granted, so the | ||
| * outcome is a pure function of the proposals and does not depend on the order the fronts are | ||
| * visited in. That is what makes the coarse grid reproducible. ---*/ | ||
| 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 is retiring and must not displace a healthy one with the | ||
| * rest of its layer. What it placed before losing stays, which is conservative. ---*/ | ||
| 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; | ||
| 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; | ||
| 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, which all-or-nothing does not allow. ---*/ | ||
| if (better(s, bids[k])) { | ||
| markFail(g); | ||
| bids[k] = s; | ||
| bidOwner[k] = f; | ||
| } else { | ||
| markFail(f); | ||
| } | ||
|
|
||
| /*--- 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; | ||
| /*--- A head-on meeting stops BOTH fronts, or the winner overshoots into the other's | ||
| * territory. A glancing contact is not a meeting and only costs the loser. ---*/ | ||
| if ((g != f) && (GeometryToolbox::DotProduct(nDim, fronts[f].dir.data(), fronts[g].dir.data()) < 0.0)) { | ||
| markFail(f); | ||
| markFail(g); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| /*--- 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); | ||
| /*--- (c) All-or-nothing acceptance: a front takes the whole layer or none of it and retires. A | ||
| * bid only becomes a claim here, so a retiring front never has to give anything back and the | ||
| * nodes it was reaching for stay available to ordinary agglomeration. ---*/ | ||
| 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] and the two carry phi. ---*/ | ||
| 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++; | ||
| 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; | ||
| if (fronts[f].failed) { | ||
| /*--- Nothing to give back: a bid only becomes a claim on acceptance. ---*/ | ||
| 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); | ||
| } | ||
|
|
||
| /*--- 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); | ||
| /*--- (d) Hand stacks across partition interfaces. A front that runs into the halo cannot go on | ||
| * here, so the footprint is sent to the owning rank, which picks the stack up and carries on. | ||
| * What crosses is the footprint, not a coarse CV, so both halves keep the same shape. It travels | ||
| * the reverse of the usual halo direction: 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; | ||
|
|
||
| Index_CoarseCV++; | ||
| any_work = true; | ||
| 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, i.e. what this rank wants the neighbour to continue. | ||
| * 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, the same way on both ranks. ---*/ | ||
| 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]; | ||
| } | ||
| } | ||
|
|
||
| pair_idx++; | ||
| if (!any_work) break; | ||
| /*--- A front that handed its WHOLE footprint over is finished here; the neighbour owns the rest | ||
| * of the stack. One that handed over only a piece keeps marching on what was left of it. ---*/ | ||
| 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); | ||
| } | ||
|
|
||
| /*--- 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; | ||
| /*--- Adopt what the neighbours sent, tags ascending so arrival order cannot change the outcome. | ||
| * A footprint whose nodes are not all free is dropped and the stack simply ends. ---*/ | ||
| 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, the same test any other layer passes. ---*/ | ||
| 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 is NOT subject to the single-layer rule the | ||
| * boundary layer gets: it opens an ordinary two-deep coarse CV and waits for its partner. ---*/ | ||
| 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; | ||
| } | ||
| if (!any_more) break; | ||
| inherited.clear(); | ||
| } |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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.
pre-commit run --allto format old commits.