From 9693757787c560b4740bdb56d4c0c714ae363f09 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Thu, 15 Oct 2026 14:18:55 +0200 Subject: [PATCH 1/8] Partition anisotropic lines onto one rank and make it the default PARMETIS_ANISOTROPY_WEIGHT (default YES) weights the ParMETIS graph edges by cell aspect ratio, gives the edges that chain points into wall-normal lines the largest weight, and moves each line to the rank holding most of its points after partitioning. Lines are then not split between ranks, which keeps the block-Jacobi ILU and line-based agglomeration from seeing fragments of lines. On stretched boundary-layer meshes this turned stalled or diverging parallel runs into converging ones (2D shear layer, 2D bump, 3D bump, flat plate with multigrid) and made iteration counts match the serial run more closely, with no case that converged without it failing with it. The load imbalance stays within 1.01-1.05 in 2D and up to about 1.3 on a large 3D mesh. The option has no effect on isotropic meshes or without ParMETIS; a warning is printed when it is set explicitly in a build without ParMETIS. Co-Authored-By: Claude Sonnet 5.5 --- Common/include/CConfig.hpp | 6 + Common/include/geometry/CPhysicalGeometry.hpp | 36 ++ Common/src/CConfig.cpp | 9 + Common/src/geometry/CPhysicalGeometry.cpp | 314 +++++++++++++++++- config_template.cfg | 5 + 5 files changed, 366 insertions(+), 4 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 9d4920c871b0..f45890acf66d 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1116,6 +1116,7 @@ class CConfig { su2double ParMETIS_tolerance; /*!< \brief Load balancing tolerance for ParMETIS. */ long ParMETIS_pointWgt; /*!< \brief Load balancing weight given to points. */ long ParMETIS_edgeWgt; /*!< \brief Load balancing weight given to edges. */ + bool ParMETIS_anisotropyWgt; /*!< \brief Weight edges by cell anisotropy, keep lines on one rank. */ unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */ bool DiscreteAdjoint, /*!< \brief AD-based discrete adjoint mode. */ DiscreteAdjointDebug; /*!< \brief Discrete adjoint debug mode using tags. */ @@ -10248,6 +10249,11 @@ class CConfig { */ long GetParMETIS_EdgeWeight() const { return ParMETIS_edgeWgt; } + /*! + * \brief Whether ParMETIS weights graph edges by cell anisotropy and keeps anisotropic lines on one rank. + */ + bool GetParMETIS_AnisotropyWeight() const { return ParMETIS_anisotropyWgt; } + /*! * \brief Find the marker index (if any) that is part of a given interface pair. * \param[in] iInterface - Number of the interface pair being tested, starting at 0. diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp index 52d8ea6da08e..e69b1a04a9b0 100644 --- a/Common/include/geometry/CPhysicalGeometry.hpp +++ b/Common/include/geometry/CPhysicalGeometry.hpp @@ -337,6 +337,42 @@ class CPhysicalGeometry final : public CGeometry { */ void PrepareAdjacency(const CConfig* config); +#if defined(HAVE_MPI) && defined(HAVE_PARMETIS) + /*! + * \brief ParMETIS edge weights that make cutting across the short edges of anisotropic cells expensive. + * \param[in] maxWeight - Cap on the weight of any edge. + * \return One weight per entry of the adjacency array. + */ + vector ComputeAnisotropyEdgeWeights(long maxWeight) const; + + /*! + * \brief Reads, for any global points of the linear partition, a value each rank holds for its own points. + * \param[in] wanted - Global indices to read. + * \param[in] localValue - Value of every local point. + * \return The value of each wanted point. + */ + vector FetchFromOwners(const vector& wanted, + const vector& localValue) const; + + /*! + * \brief Joins points into lines along the strong anisotropic edges and names each line by its lowest global index. + * \param[in] adjwgt - Anisotropy weight of every entry of the adjacency array. + * \param[out] lineEdge - Whether each entry of the adjacency array is a line edge. + * \return Global index of the root of the line of every local point, the point itself when it is on no line. + */ + vector LinkAnisotropicLines(const vector& adjwgt, vector& lineEdge) const; + + /*! + * \brief Moves each line to the partition that holds most of its points. + * \param[in] root - Line roots from LinkAnisotropicLines. + * \param[in,out] part - Partition of every local point. + * \param[out] longestLine - Number of points of the longest line. + * \return Number of points that changed partition, summed over all ranks. + */ + unsigned long KeepLinesTogether(const vector& root, vector& part, + unsigned long& longestLine) const; +#endif + /*! * \brief Find repeated nodes between two elements to identify the common face. * \param[in] first_elem - Identification of the first element. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index b0c273939be0..067331fec99e 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3095,6 +3095,9 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: ParMETIS load balancing weight for edges (equiv. to neighbors) */ addLongOption("PARMETIS_EDGE_WEIGHT", ParMETIS_edgeWgt, 1); + /* DESCRIPTION: Weight ParMETIS graph edges by cell anisotropy and keep anisotropic lines on one rank */ + addBoolOption("PARMETIS_ANISOTROPY_WEIGHT", ParMETIS_anisotropyWgt, true); + /*--- options that are used in the Hybrid RANS/LES Simulations ---*/ /*!\par CONFIG_CATEGORY:Hybrid_RANSLES Options\ingroup Config*/ @@ -4786,6 +4789,12 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i rampMUSCLValue = 1.0; } +#if !(defined(HAVE_MPI) && defined(HAVE_PARMETIS)) + if (OptionIsSet("PARMETIS_ANISOTROPY_WEIGHT") && ParMETIS_anisotropyWgt && (rank == MASTER_NODE)) { + cout << "WARNING: PARMETIS_ANISOTROPY_WEIGHT has no effect, SU2 was built without MPI and ParMETIS." << endl; + } +#endif + /*--- Check on extra Relaxation factor for Giles---*/ if(extrarelfac[1] > 0.5){ extrarelfac[1] = 0.5; diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 6928c2d922b5..3d01a1652528 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -7361,6 +7361,277 @@ void CPhysicalGeometry::SetBoundTecPlot(char mesh_filename[MAX_STRING_SIZE], boo Tecplot_File.close(); } +#if defined(HAVE_MPI) && defined(HAVE_PARMETIS) +vector CPhysicalGeometry::FetchFromOwners(const vector& wanted, + const vector& localValue) const { + CLinearPartitioner pointPartitioner(Global_nPointDomain, 0); + const unsigned long firstIndex = pointPartitioner.GetFirstIndexOnRank(rank); + + vector nSend(size, 0), nRecv(size), sendDisp(size + 1, 0), recvDisp(size + 1, 0); + for (const auto iGlobal : wanted) nSend[pointPartitioner.GetRankContainingIndex(iGlobal)]++; + SU2_MPI::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); + for (int iRank = 0; iRank < size; iRank++) { + sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; + recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; + } + + vector sendIndex(sendDisp[size]), slotOf(wanted.size()), recvIndex(recvDisp[size]); + vector next(sendDisp.begin(), sendDisp.end() - 1); + for (auto k = 0ul; k < wanted.size(); k++) { + const auto slot = next[pointPartitioner.GetRankContainingIndex(wanted[k])]++; + sendIndex[slot] = wanted[k]; + slotOf[k] = slot; + } + SU2_MPI::Alltoallv(sendIndex.data(), nSend.data(), sendDisp.data(), MPI_UNSIGNED_LONG, recvIndex.data(), nRecv.data(), + recvDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + vector reply(recvIndex.size()), answer(sendIndex.size()), result(wanted.size()); + for (auto k = 0ul; k < recvIndex.size(); k++) reply[k] = localValue[recvIndex[k] - firstIndex]; + SU2_MPI::Alltoallv(reply.data(), nRecv.data(), recvDisp.data(), MPI_UNSIGNED_LONG, answer.data(), nSend.data(), + sendDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + for (auto k = 0ul; k < wanted.size(); k++) result[k] = answer[slotOf[k]]; + return result; +} + +vector CPhysicalGeometry::LinkAnisotropicLines(const vector& adjwgt, + vector& lineEdge) const { + /*--- An edge at least this many times shorter than its element can join its points into a line. ---*/ + constexpr idx_t LINE_WEIGHT = 4; + + const unsigned long firstIndex = CLinearPartitioner(Global_nPointDomain, 0).GetFirstIndexOnRank(rank); + const vector neighbor(adjacency.begin(), adjacency.end()); + + vector strongest(nPoint, 0); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) + for (auto k = xadj[iPoint]; k < xadj[iPoint + 1]; k++) + strongest[iPoint] = max(strongest[iPoint], static_cast(adjwgt[k])); + const auto neighborStrongest = FetchFromOwners(neighbor, strongest); + + /*--- An edge joins a line if it is the strongest at one end, or nearly so at both. Equal or + * slowly varying spacing then keeps a line whole instead of splitting it where weights tie. ---*/ + lineEdge.assign(adjacency.size(), 0); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto k = xadj[iPoint]; k < xadj[iPoint + 1]; k++) { + const auto w = static_cast(adjwgt[k]); + if (w < static_cast(LINE_WEIGHT)) continue; + const bool strongestAtEnd = (w == strongest[iPoint]) || (w == neighborStrongest[k]); + const bool strongAtBoth = (2 * w >= strongest[iPoint]) && (2 * w >= neighborStrongest[k]); + lineEdge[k] = strongestAtEnd || strongAtBoth; + } + } + + /*--- Each line is named by its lowest global index, found by propagating labels along the line + * edges and pointer jumping. ---*/ + vector root(nPoint); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) root[iPoint] = firstIndex + iPoint; + for (;;) { + unsigned long changed = 0, changedGlobal = 0; + const auto neighborRoot = FetchFromOwners(neighbor, root); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto k = xadj[iPoint]; k < xadj[iPoint + 1]; k++) { + if (lineEdge[k] && (neighborRoot[k] < root[iPoint])) { + root[iPoint] = neighborRoot[k]; + changed++; + } + } + } + const auto rootOfRoot = FetchFromOwners(root, root); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + if (rootOfRoot[iPoint] < root[iPoint]) { + root[iPoint] = rootOfRoot[iPoint]; + changed++; + } + } + SU2_MPI::Allreduce(&changed, &changedGlobal, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + if (changedGlobal == 0) break; + } + return root; +} + +unsigned long CPhysicalGeometry::KeepLinesTogether(const vector& root, vector& part, + unsigned long& longestLine) const { + CLinearPartitioner pointPartitioner(Global_nPointDomain, 0); + const unsigned long firstIndex = pointPartitioner.GetFirstIndexOnRank(rank); + + /*--- The owner of each root counts in which partitions its line lies. ---*/ + vector nSend(size, 0), nRecv(size), sendDisp(size + 1, 0), recvDisp(size + 1, 0); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) nSend[pointPartitioner.GetRankContainingIndex(root[iPoint])] += 2; + SU2_MPI::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); + for (int iRank = 0; iRank < size; iRank++) { + sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; + recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; + } + vector sendBuf(sendDisp[size]), recvBuf(recvDisp[size]); + vector next(sendDisp.begin(), sendDisp.end() - 1); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + auto& slot = next[pointPartitioner.GetRankContainingIndex(root[iPoint])]; + sendBuf[slot++] = root[iPoint]; + sendBuf[slot++] = static_cast(part[iPoint]); + } + SU2_MPI::Alltoallv(sendBuf.data(), nSend.data(), sendDisp.data(), MPI_UNSIGNED_LONG, recvBuf.data(), nRecv.data(), + recvDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + map> votes; + for (auto k = 0ul; k < recvBuf.size(); k += 2) votes[recvBuf[k]][recvBuf[k + 1]]++; + + unsigned long longest = 0; + for (const auto& line : votes) { + unsigned long length = 0; + for (const auto& vote : line.second) length += vote.second; + longest = max(longest, length); + } + SU2_MPI::Allreduce(&longest, &longestLine, 1, MPI_UNSIGNED_LONG, MPI_MAX, SU2_MPI::GetComm()); + + vector decision(nPoint); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) decision[iPoint] = static_cast(part[iPoint]); + for (const auto& line : votes) { + unsigned long best = 0, bestCount = 0; + for (const auto& vote : line.second) { + if (vote.second > bestCount) { + best = vote.first; + bestCount = vote.second; + } + } + decision[line.first - firstIndex] = best; + } + + /*--- Every point takes the partition of its line, unless that would leave a rank without points. ---*/ + const auto chosen = FetchFromOwners(root, decision); + vector count(size, 0), countGlobal(size); + for (const auto iPart : chosen) count[iPart]++; + SU2_MPI::Allreduce(count.data(), countGlobal.data(), size, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + if (*min_element(countGlobal.begin(), countGlobal.end()) == 0) return 0; + + unsigned long moved = 0, movedGlobal = 0; + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + if (static_cast(part[iPoint]) != chosen[iPoint]) { + part[iPoint] = static_cast(chosen[iPoint]); + moved++; + } + } + SU2_MPI::Allreduce(&moved, &movedGlobal, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + return movedGlobal; +} + +/*! + * \brief Compute graph edge weights that encode mesh anisotropy for partitioning. + * + * The routine builds edge weights from geometric information of neighboring points. + * Since element connectivity can reference points owned by other MPI ranks, it first + * exchanges the remote coordinates that are required locally. The resulting weights are + * scaled/clamped to the user-provided \p maxWeight and returned in graph-edge order. + * + * \param[in] maxWeight Upper bound used when converting anisotropy metrics to integer + * partitioning weights. + * \return Edge weights suitable for graph partitioning (one value per graph edge entry). + */ +vector CPhysicalGeometry::ComputeAnisotropyEdgeWeights(long maxWeight) const { + /*--- Build local/global point-index helpers used throughout communication. ---*/ + CLinearPartitioner pointPartitioner(Global_nPointDomain, 0); + const unsigned long firstIndex = pointPartitioner.GetFirstIndexOnRank(rank); + const auto isLocal = [&](unsigned long iGlobal) { + return (iGlobal >= firstIndex) && (iGlobal - firstIndex < nPoint); + }; + + /*--- Identify remote element nodes whose coordinates must be requested. ---*/ + vector> request(size); + for (auto iElem = 0ul; iElem < nElem; iElem++) { + for (auto iNode = 0u; iNode < elem[iElem]->GetnNodes(); iNode++) { + const auto iGlobal = elem[iElem]->GetNode(iNode); + if (!isLocal(iGlobal)) request[pointPartitioner.GetRankContainingIndex(iGlobal)].push_back(iGlobal); + } + } + + /*--- Exchange request sizes (after deduplication) to prepare all-to-all transfers. ---*/ + vector nSend(size), nRecv(size), sendDisp(size + 1, 0), recvDisp(size + 1, 0); + for (int iRank = 0; iRank < size; iRank++) { + auto& list = request[iRank]; + sort(list.begin(), list.end()); + list.erase(unique(list.begin(), list.end()), list.end()); + nSend[iRank] = static_cast(list.size()); + } + SU2_MPI::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); + + for (int iRank = 0; iRank < size; iRank++) { + sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; + recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; + } + + vector sendIndex, recvIndex(recvDisp[size]); + sendIndex.reserve(sendDisp[size]); + for (const auto& list : request) sendIndex.insert(sendIndex.end(), list.begin(), list.end()); + + SU2_MPI::Alltoallv(sendIndex.data(), nSend.data(), sendDisp.data(), MPI_UNSIGNED_LONG, recvIndex.data(), nRecv.data(), + recvDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + vector replyCoord(recvIndex.size() * nDim), remoteCoord(sendIndex.size() * nDim); + for (auto i = 0ul; i < recvIndex.size(); i++) + for (auto iDim = 0u; iDim < nDim; iDim++) + replyCoord[i * nDim + iDim] = SU2_TYPE::GetValue(nodes->GetCoord(recvIndex[i] - firstIndex, iDim)); + + for (int iRank = 0; iRank < size; iRank++) { + nSend[iRank] *= nDim; + nRecv[iRank] *= nDim; + sendDisp[iRank] *= nDim; + recvDisp[iRank] *= nDim; + } + SU2_MPI::Alltoallv(replyCoord.data(), nRecv.data(), recvDisp.data(), MPI_DOUBLE, remoteCoord.data(), nSend.data(), + sendDisp.data(), MPI_DOUBLE, SU2_MPI::GetComm()); + + unordered_map remoteSlot; + for (auto i = 0ul; i < sendIndex.size(); i++) remoteSlot[sendIndex[i]] = i; + + auto getCoord = [&](unsigned long iGlobal, unsigned short iDim) -> passivedouble { + if (isLocal(iGlobal)) return SU2_TYPE::GetValue(nodes->GetCoord(iGlobal - firstIndex, iDim)); + return remoteCoord[remoteSlot.at(iGlobal) * nDim + iDim]; + }; + + /*--- Weight is the longest node distance of the element over the edge length, maximized over elements. ---*/ + + vector weight(adjacency.size(), 1); + + for (auto iElem = 0ul; iElem < nElem; iElem++) { + const auto nNodes = elem[iElem]->GetnNodes(); + passivedouble coord[N_POINTS_HEXAHEDRON][MAXNDIM] = {}; + for (auto iNode = 0u; iNode < nNodes; iNode++) + for (auto iDim = 0u; iDim < nDim; iDim++) coord[iNode][iDim] = getCoord(elem[iElem]->GetNode(iNode), iDim); + + auto distance2 = [&](unsigned short a, unsigned short b) { + passivedouble d2 = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) d2 += pow(coord[a][iDim] - coord[b][iDim], 2); + return d2; + }; + + passivedouble maxDist2 = 0.0; + for (auto a = 0u; a < nNodes; a++) + for (auto b = a + 1; b < nNodes; b++) maxDist2 = max(maxDist2, distance2(a, b)); + + for (auto a = 0u; a < nNodes; a++) { + const auto iGlobal = elem[iElem]->GetNode(a); + if (!isLocal(iGlobal)) continue; + const auto iLocal = iGlobal - firstIndex; + const auto begin = adjacency.begin() + xadj[iLocal]; + const auto end = adjacency.begin() + xadj[iLocal + 1]; + + for (auto b = 0u; b < nNodes; b++) { + const auto d2 = distance2(a, b); + if ((b == a) || (d2 <= 0.0)) continue; + + const auto jGlobal = static_cast(elem[iElem]->GetNode(b)); + const auto it = lower_bound(begin, end, jGlobal); + if ((it == end) || (*it != jGlobal)) continue; + + const auto ratio = static_cast(sqrt(maxDist2 / d2)); + auto& w = weight[it - adjacency.begin()]; + w = max(w, min(maxWeight, max(1, ratio))); + } + } + } + return weight; +} +#endif + void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { /*--- We need to have parallel support with MPI and have the ParMETIS library compiled and linked for parallel graph partitioning. ---*/ @@ -7412,6 +7683,23 @@ void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { vwgt[iPoint] = wp + we * (xadj[iPoint + 1] - xadj[iPoint]); } + /*--- Keep partition boundaries off the short edges of anisotropic cells. The edges that chain + * points into wall-normal lines get the largest weight, so the graph coarsening merges each + * line before it balances. The cap keeps the total edge weight within the range of idx_t. ---*/ + + const bool anisotropic = config->GetParMETIS_AnisotropyWeight(); + vector adjwgt; + vector lineRoot; + if (anisotropic) { + constexpr long maxWeight = (sizeof(idx_t) >= 8) ? 1000000 : 1000; + adjwgt = ComputeAnisotropyEdgeWeights(maxWeight); + vector lineEdge; + lineRoot = LinkAnisotropicLines(adjwgt, lineEdge); + for (auto k = 0ul; k < adjwgt.size(); k++) + if (lineEdge[k]) adjwgt[k] = maxWeight; + wgtflag = 3; + } + /*--- Create some structures that ParMETIS needs to output the partitioning. ---*/ idx_t edgecut; @@ -7420,12 +7708,30 @@ void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { /*--- Calling ParMETIS ---*/ if (rank == MASTER_NODE) cout << "Calling ParMETIS..."; - auto err = - ParMETIS_V3_PartKway(vtxdist.data(), xadj.data(), adjacency.data(), vwgt.data(), nullptr, &wgtflag, &numflag, - &ncon, &nparts, tpwgts.data(), &ubvec, options, &edgecut, part.data(), &comm); + const int err = ParMETIS_V3_PartKway(vtxdist.data(), xadj.data(), adjacency.data(), vwgt.data(), + anisotropic ? adjwgt.data() : nullptr, &wgtflag, &numflag, &ncon, &nparts, + tpwgts.data(), &ubvec, options, &edgecut, part.data(), &comm); if (err != METIS_OK) SU2_MPI::Error("Partitioning failed.", CURRENT_FUNCTION); if (rank == MASTER_NODE) { - cout << " graph partitioning complete (" << edgecut << " edge cuts)." << endl; + cout << " graph partitioning complete (" << edgecut << (anisotropic ? " weighted" : "") << " edge cuts)." << endl; + } + + /*--- A line split between ranks loses its strongest couplings in the preconditioner, and its + * fragments give tangled coarse control volumes. ---*/ + + if (anisotropic) { + unsigned long longestLine = 0; + const auto moved = KeepLinesTogether(lineRoot, part, longestLine); + /*--- The imbalance is measured in the vertex weights ParMETIS balances. ---*/ + vector load(size, 0), loadGlobal(size); + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) load[part[iPoint]] += vwgt[iPoint]; + SU2_MPI::Allreduce(load.data(), loadGlobal.data(), size, MPI_UNSIGNED_LONG, MPI_SUM, comm); + const auto maxLoad = *max_element(loadGlobal.begin(), loadGlobal.end()); + const auto sumLoad = accumulate(loadGlobal.begin(), loadGlobal.end(), 0ul); + if (rank == MASTER_NODE) { + cout << "Longest anisotropic line " << longestLine << " points, moved " << moved + << " points to keep lines on one rank, load imbalance " << double(maxLoad) * size / sumLoad << "." << endl; + } } /*--- Store the results of the partitioning (note that this is local diff --git a/config_template.cfg b/config_template.cfg index 5e611f22ed48..33ea25d39891 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -2518,6 +2518,11 @@ PARMETIS_TOLERANCE= 0.02 PARMETIS_EDGE_WEIGHT= 1 PARMETIS_POINT_WEIGHT= 0 % +% Weight the partitioning graph by cell anisotropy and keep chains of strongly coupled points +% (wall-normal lines in stretched boundary layers) on a single rank, so that line-based +% agglomeration and the block-Jacobi preconditioner do not see split lines. +PARMETIS_ANISOTROPY_WEIGHT= YES +% % ----------------------- SOBOLEV GRADIENT SMOOTHING OPTIONS ----------------------% % % Activate the gradient smoothing solver for the discrete adjoint driver (NO, YES) From 4cadbd9c9b9020212969c50533ce0038d5de08aa Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 17 Oct 2026 18:22:49 +0200 Subject: [PATCH 2/8] precommit --- .pre-commit-config.yaml | 2 +- Common/src/geometry/CPhysicalGeometry.cpp | 22 +++++------ SU2_PY/FSI_tools/FSIInterface.py | 12 +++--- SU2_PY/SU2/eval/functions.py | 1 - SU2_PY/SU2/io/filelock.py | 1 + SU2_PY/SU2/io/redirect.py | 1 + SU2_PY/SU2/io/state.py | 1 - SU2_PY/SU2/io/tools.py | 10 +++-- SU2_PY/SU2/opt/scipy_tools.py | 1 - SU2_PY/SU2/run/deform.py | 1 - SU2_PY/SU2/run/direct.py | 6 +-- SU2_PY/SU2/run/projection.py | 1 - SU2_PY/SU2/util/bunch.py | 48 +++++++++++------------ SU2_PY/SU2/util/filter_adjoint.py | 2 +- SU2_PY/SU2/util/ordered_bunch.py | 48 +++++++++++------------ SU2_PY/SU2/util/ordered_dict.py | 6 +-- SU2_PY/aeroacoustics.py | 1 - SU2_PY/tests/test_aeroacoustics.py | 1 - SU2_PY/topology_optimization.py | 1 + meson.py | 1 - 20 files changed, 83 insertions(+), 84 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 1e55bb83a103..a9e9f8d307ce 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -18,7 +18,7 @@ repos: types_or: [c++, c] # black repo for python formatting - repo: https://github.com/ambv/black - rev: 22.6.0 + rev: 26.10.1 hooks: - id: black exclude: | diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 3d01a1652528..ef7de89ec30f 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -2426,9 +2426,9 @@ void CPhysicalGeometry::LoadSurfaceElements(CConfig* config, CGeometry* geometry unsigned long Local_Nodes[N_POINTS_HEXAHEDRON]; - vector > Line_List; - vector > BoundTria_List; - vector > BoundQuad_List; + vector> Line_List; + vector> BoundTria_List; + vector> BoundQuad_List; vector Marker_Local; @@ -2967,14 +2967,14 @@ void CPhysicalGeometry::SetSendReceive(const CConfig* config) { unsigned short iNode, jNode; vector::iterator it; - vector > + vector> SendTransfLocal; /*!< \brief Vector to store the type of transformation for this send point. */ - vector > + vector> ReceivedTransfLocal; /*!< \brief Vector to store the type of transformation for this received point. */ - vector > SendDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the point - index of the node that must me sended. */ - vector > ReceivedDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the - point index of the node that must me sended. */ + vector> SendDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the point + index of the node that must me sended. */ + vector> ReceivedDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the + point index of the node that must me sended. */ unordered_map::const_iterator MI; @@ -4426,7 +4426,7 @@ void CPhysicalGeometry::SetPositive_ZArea(CConfig* config) { } void CPhysicalGeometry::SetPoint_Connectivity() { - vector > points(nPoint); + vector> points(nPoint); SU2_OMP_PARALLEL { unsigned short Node_Neighbor, iNode, iNeighbor; @@ -4434,7 +4434,7 @@ void CPhysicalGeometry::SetPoint_Connectivity() { /*--- Loop over all the elements ---*/ BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - vector > elems(nPoint); + vector> elems(nPoint); for (iElem = 0; iElem < nElem; iElem++) { /*--- Loop over all the nodes of an element ---*/ diff --git a/SU2_PY/FSI_tools/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py index 01f2fa2e14bc..bf5eabe16dc2 100644 --- a/SU2_PY/FSI_tools/FSIInterface.py +++ b/SU2_PY/FSI_tools/FSIInterface.py @@ -478,9 +478,9 @@ def connect(self, FSI_config, FluidSolver, SolidSolver): ) del sendBuffPhysical else: - self.fluidPhysicalInterfaceNodesDistribution[ - 0 - ] = self.nFluidInterfacePhysicalNodes + self.fluidPhysicalInterfaceNodesDistribution[0] = ( + self.nFluidInterfacePhysicalNodes + ) # Same thing for the solid part self.solidPhysicalInterfaceNodesDistribution = np.zeros(MPIsize, dtype=int) @@ -491,9 +491,9 @@ def connect(self, FSI_config, FluidSolver, SolidSolver): ) del sendBuffPhysical else: - self.solidPhysicalInterfaceNodesDistribution[ - 0 - ] = self.nSolidInterfacePhysicalNodes + self.solidPhysicalInterfaceNodesDistribution[0] = ( + self.nSolidInterfacePhysicalNodes + ) # --- Calculate and store the global indexing of interface physical nodes on each processor and allgather the information --- if self.have_MPI: diff --git a/SU2_PY/SU2/eval/functions.py b/SU2_PY/SU2/eval/functions.py index fd6c1103a519..7f6bf4ee41cf 100644 --- a/SU2_PY/SU2/eval/functions.py +++ b/SU2_PY/SU2/eval/functions.py @@ -35,7 +35,6 @@ from .. import util as su2util from ..io import redirect_folder, redirect_output - # ---------------------------------------------------------------------- # Main Function Interface # ---------------------------------------------------------------------- diff --git a/SU2_PY/SU2/io/filelock.py b/SU2_PY/SU2/io/filelock.py index 45190e1711a0..e1101031e78a 100644 --- a/SU2_PY/SU2/io/filelock.py +++ b/SU2_PY/SU2/io/filelock.py @@ -28,6 +28,7 @@ import os, time, errno from random import random + # ------------------------------------------------------------------- # File Lock Class # ------------------------------------------------------------------- diff --git a/SU2_PY/SU2/io/redirect.py b/SU2_PY/SU2/io/redirect.py index 161d4cc31f0c..ff1b7a5d4ea1 100644 --- a/SU2_PY/SU2/io/redirect.py +++ b/SU2_PY/SU2/io/redirect.py @@ -32,6 +32,7 @@ import os, sys, shutil, copy, glob from .tools import add_suffix, make_link, expand_part + # ------------------------------------------------------------------- # Output Redirection # ------------------------------------------------------------------- diff --git a/SU2_PY/SU2/io/state.py b/SU2_PY/SU2/io/state.py index f67a9004b81d..d904258e8e60 100644 --- a/SU2_PY/SU2/io/state.py +++ b/SU2_PY/SU2/io/state.py @@ -44,7 +44,6 @@ from ..util import bunch from ..util import ordered_bunch - # ---------------------------------------------------------------------- # State Factory # ---------------------------------------------------------------------- diff --git a/SU2_PY/SU2/io/tools.py b/SU2_PY/SU2/io/tools.py index aaeccf9433fc..237ac021002f 100755 --- a/SU2_PY/SU2/io/tools.py +++ b/SU2_PY/SU2/io/tools.py @@ -302,6 +302,7 @@ def getTurboPerfIndex(nZones=1): "EFFICIENCY": "CL/CD", } + # ------------------------------------------------------------------- # Include per-surface output from History File # ------------------------------------------------------------------- @@ -361,9 +362,9 @@ def read_aerodynamics( historyOutFields[this_objfun]["TYPE"] == "COEFFICIENT" or historyOutFields[this_objfun]["TYPE"] == "D_COEFFICIENT" ): - Func_Values[ - this_objfun + "[" + str(iZone) + "]" - ] = history_data[this_objfun + "[" + str(iZone) + "]"] + Func_Values[this_objfun + "[" + str(iZone) + "]"] = ( + history_data[this_objfun + "[" + str(iZone) + "]"] + ) if "TIME_MARCHING" in special_cases: # for unsteady cases, average time-accurate objective function values @@ -626,6 +627,7 @@ def get_dvMap(): #: def get_dvMap() + # ------------------------------------------------------------------- # Get Design Variable Kind Name from ID # ------------------------------------------------------------------- @@ -640,6 +642,7 @@ def get_dvKind(kindID): # def get_dvKind() + # ------------------------------------------------------------------- # Get Design Variable Kind ID from Name # ------------------------------------------------------------------- @@ -996,6 +999,7 @@ def get_specialCases(config): #: def get_specialCases() + # ------------------------------------------------------------------- # Check Fluid Structure Interaction # ------------------------------------------------------------------- diff --git a/SU2_PY/SU2/opt/scipy_tools.py b/SU2_PY/SU2/opt/scipy_tools.py index e8551b67ac08..1029bc6f7df4 100644 --- a/SU2_PY/SU2/opt/scipy_tools.py +++ b/SU2_PY/SU2/opt/scipy_tools.py @@ -34,7 +34,6 @@ from .. import eval as su2eval from numpy import array, zeros - # ------------------------------------------------------------------- # Scipy SLSQP # ------------------------------------------------------------------- diff --git a/SU2_PY/SU2/run/deform.py b/SU2_PY/SU2/run/deform.py index 8ead56f11509..7feb9feaef91 100644 --- a/SU2_PY/SU2/run/deform.py +++ b/SU2_PY/SU2/run/deform.py @@ -34,7 +34,6 @@ from .. import io as su2io from .interface import DEF as SU2_DEF - # ---------------------------------------------------------------------- # Mesh Deformation # ---------------------------------------------------------------------- diff --git a/SU2_PY/SU2/run/direct.py b/SU2_PY/SU2/run/direct.py index 60825cae7a63..f6f91089e513 100644 --- a/SU2_PY/SU2/run/direct.py +++ b/SU2_PY/SU2/run/direct.py @@ -93,9 +93,9 @@ def direct(config): # check for 'RESTART_ITER' is to avoid forced restart situation in "compute_polar.py"... if konfig.get("RESTART_SOL", "NO") == "YES" and konfig.get("RESTART_ITER", 1) != 1: if konfig.get("CONFIG_LIST", []) != []: - konfig[ - "CONV_FILENAME" - ] = "config_CFD" # master cfg is always config_CFD. Hardcoded names are prob nt ideal. + konfig["CONV_FILENAME"] = ( + "config_CFD" # master cfg is always config_CFD. Hardcoded names are prob nt ideal. + ) restart_iter = "_" + str(konfig["RESTART_ITER"]).zfill(5) history_filename = konfig["CONV_FILENAME"] + restart_iter + plot_extension else: diff --git a/SU2_PY/SU2/run/projection.py b/SU2_PY/SU2/run/projection.py index 80741ba70759..4811c947a6a1 100644 --- a/SU2_PY/SU2/run/projection.py +++ b/SU2_PY/SU2/run/projection.py @@ -35,7 +35,6 @@ from .. import util as su2util from .interface import DOT as SU2_DOT - # ---------------------------------------------------------------------- # Gradient Projection # ---------------------------------------------------------------------- diff --git a/SU2_PY/SU2/util/bunch.py b/SU2_PY/SU2/util/bunch.py index 156e681e7f78..44afe34f8eee 100644 --- a/SU2_PY/SU2/util/bunch.py +++ b/SU2_PY/SU2/util/bunch.py @@ -1,29 +1,29 @@ #!/usr/bin/env python -""" Bunch is a subclass of dict with attribute-style access. - - >>> b = Bunch() - >>> b.hello = 'world' - >>> b.hello - 'world' - >>> b['hello'] += "!" - >>> b.hello - 'world!' - >>> b.foo = Bunch(lol=True) - >>> b.foo.lol - True - >>> b.foo is b['foo'] - True - - It is safe to import * from this module: - - __all__ = ('Bunch', 'bunchify','unbunchify') - - un/bunchify provide dictionary conversion; Bunches can also be - converted via Bunch.to/fromDict(). - - original source: - https://pypi.python.org/pypi/bunch +"""Bunch is a subclass of dict with attribute-style access. + +>>> b = Bunch() +>>> b.hello = 'world' +>>> b.hello +'world' +>>> b['hello'] += "!" +>>> b.hello +'world!' +>>> b.foo = Bunch(lol=True) +>>> b.foo.lol +True +>>> b.foo is b['foo'] +True + +It is safe to import * from this module: + + __all__ = ('Bunch', 'bunchify','unbunchify') + +un/bunchify provide dictionary conversion; Bunches can also be +converted via Bunch.to/fromDict(). + +original source: +https://pypi.python.org/pypi/bunch """ diff --git a/SU2_PY/SU2/util/filter_adjoint.py b/SU2_PY/SU2/util/filter_adjoint.py index d3a541bbc6bf..c290812786bf 100644 --- a/SU2_PY/SU2/util/filter_adjoint.py +++ b/SU2_PY/SU2/util/filter_adjoint.py @@ -79,7 +79,7 @@ def main(): metavar="CHORD", ) - (options, args) = parser.parse_args() + options, args = parser.parse_args() options.chord_length = float(options.chord_length) # run filter diff --git a/SU2_PY/SU2/util/ordered_bunch.py b/SU2_PY/SU2/util/ordered_bunch.py index 7aa0f33427de..df8c4bf26172 100644 --- a/SU2_PY/SU2/util/ordered_bunch.py +++ b/SU2_PY/SU2/util/ordered_bunch.py @@ -1,29 +1,29 @@ #!/usr/bin/env python -""" OrderedBunch is a subclass of OrderedDict with attribute-style access. - - >>> b = OrderedBunch() - >>> b.hello = 'world' - >>> b.hello - 'world' - >>> b['hello'] += "!" - >>> b.hello - 'world!' - >>> b.foo = OrderedBunch(lol=True) - >>> b.foo.lol - True - >>> b.foo is b['foo'] - True - - It is safe to import * from this module: - - __all__ = ('OrderedBunch', 'ordered_bunchify','ordered_unbunchify') - - ordered_un/bunchify provide dictionary conversion; Bunches can also be - converted via OrderedBunch.to/fromOrderedDict(). - - original source: - https://pypi.python.org/pypi/bunch +"""OrderedBunch is a subclass of OrderedDict with attribute-style access. + +>>> b = OrderedBunch() +>>> b.hello = 'world' +>>> b.hello +'world' +>>> b['hello'] += "!" +>>> b.hello +'world!' +>>> b.foo = OrderedBunch(lol=True) +>>> b.foo.lol +True +>>> b.foo is b['foo'] +True + +It is safe to import * from this module: + + __all__ = ('OrderedBunch', 'ordered_bunchify','ordered_unbunchify') + +ordered_un/bunchify provide dictionary conversion; Bunches can also be +converted via OrderedBunch.to/fromOrderedDict(). + +original source: +https://pypi.python.org/pypi/bunch """ from .ordered_dict import OrderedDict diff --git a/SU2_PY/SU2/util/ordered_dict.py b/SU2_PY/SU2/util/ordered_dict.py index 776aa1ca55fc..b4fbe21066f7 100644 --- a/SU2_PY/SU2/util/ordered_dict.py +++ b/SU2_PY/SU2/util/ordered_dict.py @@ -1,6 +1,6 @@ -""" Backport of OrderedDict() class that runs on Python 2.4, 2.5, 2.6, 2.7 and pypy. - Passes Python2.7's test suite and incorporates all the latest updates. - {{{ http://code.activestate.com/recipes/576693/ (r9) +"""Backport of OrderedDict() class that runs on Python 2.4, 2.5, 2.6, 2.7 and pypy. +Passes Python2.7's test suite and incorporates all the latest updates. +{{{ http://code.activestate.com/recipes/576693/ (r9) """ try: diff --git a/SU2_PY/aeroacoustics.py b/SU2_PY/aeroacoustics.py index 5243a6668f6e..ae2aca0f4bba 100755 --- a/SU2_PY/aeroacoustics.py +++ b/SU2_PY/aeroacoustics.py @@ -35,7 +35,6 @@ import numpy as np - WINDOWS = { "hann": np.hanning, "hamming": np.hamming, diff --git a/SU2_PY/tests/test_aeroacoustics.py b/SU2_PY/tests/test_aeroacoustics.py index 31cc242aba7c..573141fcf496 100644 --- a/SU2_PY/tests/test_aeroacoustics.py +++ b/SU2_PY/tests/test_aeroacoustics.py @@ -10,7 +10,6 @@ import numpy as np - MODULE_PATH = Path(__file__).parents[1] / "aeroacoustics.py" SPEC = importlib.util.spec_from_file_location("aeroacoustics", MODULE_PATH) AEROACOUSTICS = importlib.util.module_from_spec(SPEC) diff --git a/SU2_PY/topology_optimization.py b/SU2_PY/topology_optimization.py index 068e006cb878..b8d57c1cebd5 100755 --- a/SU2_PY/topology_optimization.py +++ b/SU2_PY/topology_optimization.py @@ -225,6 +225,7 @@ def con_der(self, x): ####### Helpers ####### + # updates the parameters in the config files def update_settings(fnames, params): for fname in fnames: diff --git a/meson.py b/meson.py index ca03e47b3f0c..467dce5055d5 100755 --- a/meson.py +++ b/meson.py @@ -30,7 +30,6 @@ sys.path.append(sys.path[0]) import preconfigure - if __name__ == "__main__": if sys.version_info[0] < 3: raise Exception("Script must be run using Python 3") From af9cf7f4d4fea595fe1b439337e6334cd28f0d2f Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 18 Oct 2026 09:06:26 +0200 Subject: [PATCH 3/8] Use su2double buffers for the coordinate exchange in the anisotropy weights In AD builds MPI_DOUBLE maps to AMPI_ADOUBLE, which expects active su2double buffers. The passivedouble buffers were overrun, corrupting the heap and crashing the AD regression cases that partition with ParMETIS. Co-Authored-By: Claude Sonnet 5.5 --- Common/src/geometry/CPhysicalGeometry.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index ef7de89ec30f..382e5a6e5748 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -7565,10 +7565,11 @@ vector CPhysicalGeometry::ComputeAnisotropyEdgeWeights(long maxWeight) co SU2_MPI::Alltoallv(sendIndex.data(), nSend.data(), sendDisp.data(), MPI_UNSIGNED_LONG, recvIndex.data(), nRecv.data(), recvDisp.data(), MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); - vector replyCoord(recvIndex.size() * nDim), remoteCoord(sendIndex.size() * nDim); + /*--- MPI_DOUBLE is an active type in AD builds, so the buffers hold su2double. ---*/ + vector replyCoord(recvIndex.size() * nDim), remoteCoord(sendIndex.size() * nDim); for (auto i = 0ul; i < recvIndex.size(); i++) for (auto iDim = 0u; iDim < nDim; iDim++) - replyCoord[i * nDim + iDim] = SU2_TYPE::GetValue(nodes->GetCoord(recvIndex[i] - firstIndex, iDim)); + replyCoord[i * nDim + iDim] = nodes->GetCoord(recvIndex[i] - firstIndex, iDim); for (int iRank = 0; iRank < size; iRank++) { nSend[iRank] *= nDim; @@ -7584,7 +7585,7 @@ vector CPhysicalGeometry::ComputeAnisotropyEdgeWeights(long maxWeight) co auto getCoord = [&](unsigned long iGlobal, unsigned short iDim) -> passivedouble { if (isLocal(iGlobal)) return SU2_TYPE::GetValue(nodes->GetCoord(iGlobal - firstIndex, iDim)); - return remoteCoord[remoteSlot.at(iGlobal) * nDim + iDim]; + return SU2_TYPE::GetValue(remoteCoord[remoteSlot.at(iGlobal) * nDim + iDim]); }; /*--- Weight is the longest node distance of the element over the edge length, maximized over elements. ---*/ From e687638234c9d11d72afb421aaf5d896801fd346 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 18 Oct 2026 19:33:37 +0200 Subject: [PATCH 4/8] Sum the streamwise periodic pressure drop adjoint over all ranks Every rank holds a copy of the pressure drop, so the adjoint extracted on one rank was only its local contribution and Adjoint_SolExtra depended on the partition. The extracted adjoint is now summed over the ranks, and the seed is divided by the number of ranks because the reverse Allreduce of the massflow sums the seeds of all copies. Co-Authored-By: Claude Opus 5.5 --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 11b4e3f9f063..cc832823d7b8 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -3665,14 +3665,21 @@ unsigned long CIncEulerSolver::RegisterSolutionExtra(bool input, const CConfig* void CIncEulerSolver::SetAdjoint_SolutionExtra(const su2activevector& adj_sol, const CConfig* config) { SU2_ZONE_SCOPED if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { - SU2_TYPE::SetDerivative(SPvalsUpdated.Streamwise_Periodic_PressureDrop, SU2_TYPE::GetValue(adj_sol[0])); + /*--- Every rank holds a copy of the pressure drop, the reverse Allreduce of the massflow sums their seeds. ---*/ + SU2_TYPE::SetDerivative(SPvalsUpdated.Streamwise_Periodic_PressureDrop, SU2_TYPE::GetValue(adj_sol[0]) / size); } } void CIncEulerSolver::ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, const CConfig* config) { SU2_ZONE_SCOPED if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { - adj_sol[0] = SU2_TYPE::GetDerivative(SPvals.Streamwise_Periodic_PressureDrop); + /*--- Sum the adjoints of the copies of the pressure drop on all ranks. ---*/ + SU2_OMP_SAFE_GLOBAL_ACCESS( + const passivedouble local = SU2_TYPE::GetDerivative(SPvals.Streamwise_Periodic_PressureDrop); + passivedouble global = 0.0; + SelectMPIWrapper::W::Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + adj_sol[0] = global; + ) } } From b93f3bfb16501fb3490a59c4b0ce75bff039c0e4 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 18 Oct 2026 19:33:37 +0200 Subject: [PATCH 5/8] Copy the owner rows into the halo rows before transposing the Jacobian TransposeInPlace turns the blocks of halo rows into owned rows of the transposed matrix. Boundary conditions modify only owned rows (symmetry planes project out the normal momentum, walls delete velocity, energy and turbulence rows), so wherever a partition boundary met such a boundary the adjoint linear system had spurious couplings. The primal never reads halo rows and was unaffected. On da_sp_pinArray_cht_2d_mf at CFL 100 this made the adjoint diverge on some partitions (2 ranks with anisotropic ParMETIS weights, 3 and 4 ranks without), with the adjoint linear solver reaching only 1e-1.4 instead of 1e-3. All partitions from 1 to 4 ranks now give the same results. The reference values of that case are updated, they include the pressure drop adjoint fix. Co-Authored-By: Claude Opus 5.5 --- Common/include/linear_algebra/CSysMatrix.hpp | 6 ++ Common/src/linear_algebra/CSysMatrix.cpp | 63 ++++++++++++++++++++ Common/src/linear_algebra/CSysSolve.cpp | 8 ++- TestCases/parallel_regression_AD.py | 2 +- 4 files changed, 77 insertions(+), 2 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 70dd943690ec..1adc0187b4ee 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -1293,6 +1293,12 @@ class CSysMatrix { */ void SetDiagonalAsColumnSum(); + /*! + * \brief Overwrites the halo rows with the rows of the owning ranks, which hold the boundary condition modifications. + * \param[in] geometry - Geometry with the point-to-point communication pattern. + */ + void CopyHaloRowsFromOwners(const CGeometry* geometry); + /*! * \brief Transposes the matrix, any preconditioner that was computed may be invalid. */ diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 4e389fd0b3dc..8a25e085cb41 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -1706,6 +1706,69 @@ void CSysMatrix::SetDiagonalAsColumnSum() { END_SU2_OMP_FOR } +template +void CSysMatrix::CopyHaloRowsFromOwners(const CGeometry* geometry) { + SU2_ZONE_SCOPED + const int size = SU2_MPI::GetSize(); + if (size == 1) return; + + using MPIWrapper = SelectMPIWrapper::W; + const auto* nodes = geometry->nodes; + const auto blkSize = nVar * nEqn; + + /*--- Each point sent to a rank is followed by its blocks in the columns the sender does not own, + * as the block count, then the global column index and the block for each. ---*/ + vector> sendTo(size); + for (int iMsg = 0; iMsg < geometry->nP2PSend; iMsg++) { + auto& buf = sendTo[geometry->Neighbors_P2PSend[iMsg]]; + for (auto k = geometry->nPoint_P2PSend[iMsg]; k < geometry->nPoint_P2PSend[iMsg + 1]; k++) { + const auto iPoint = geometry->Local_Point_P2PSend[k]; + const auto countPos = buf.size(); + buf.push_back(0); + for (const auto jPoint : nodes->GetPoints(iPoint)) { + if (nodes->GetDomain(jPoint)) continue; + const auto* blk = GetBlock(iPoint, jPoint); + if (!blk) continue; + buf.push_back(nodes->GetGlobalIndex(jPoint)); + for (auto i = 0ul; i < blkSize; i++) buf.push_back(SU2_TYPE::GetValue(blk[i])); + buf[countPos] += 1; + } + } + } + + vector nSend(size), nRecv(size), sendDisp(size + 1, 0), recvDisp(size + 1, 0); + for (int iRank = 0; iRank < size; iRank++) nSend[iRank] = sendTo[iRank].size(); + MPIWrapper::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); + for (int iRank = 0; iRank < size; iRank++) { + sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; + recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; + } + vector sendBuf(sendDisp[size]), recvBuf(recvDisp[size]); + for (int iRank = 0; iRank < size; iRank++) + copy(sendTo[iRank].begin(), sendTo[iRank].end(), sendBuf.begin() + sendDisp[iRank]); + MPIWrapper::Alltoallv(sendBuf.data(), nSend.data(), sendDisp.data(), MPI_DOUBLE, recvBuf.data(), nRecv.data(), + recvDisp.data(), MPI_DOUBLE, SU2_MPI::GetComm()); + + /*--- The points of a message arrive in the order of the receive list from that rank. Only the + * blocks in owned columns are copied, these become owned rows when the matrix is transposed. ---*/ + for (int iMsg = 0; iMsg < geometry->nP2PRecv; iMsg++) { + const auto* data = recvBuf.data() + recvDisp[geometry->Neighbors_P2PRecv[iMsg]]; + for (auto k = geometry->nPoint_P2PRecv[iMsg]; k < geometry->nPoint_P2PRecv[iMsg + 1]; k++) { + const auto iPoint = geometry->Local_Point_P2PRecv[k]; + const auto nBlk = static_cast(*data++); + for (auto iBlk = 0ul; iBlk < nBlk; iBlk++, data += blkSize) { + const auto jGlobal = static_cast(*data++); + for (const auto jPoint : nodes->GetPoints(iPoint)) { + if (!nodes->GetDomain(jPoint) || nodes->GetGlobalIndex(jPoint) != jGlobal) continue; + auto* blk = GetBlock(iPoint, jPoint); + if (blk) for (auto i = 0ul; i < blkSize; i++) blk[i] = data[i]; + break; + } + } + } + } +} + template void CSysMatrix::TransposeInPlace() { SU2_ZONE_SCOPED diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index c1a356dfbcce..2ed0205f9a0b 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -1590,7 +1590,11 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con /*--- Build preconditioner for the transposed Jacobian ---*/ - if (RequiresTranspose) Jacobian.TransposeInPlace(); + if (RequiresTranspose) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { Jacobian.CopyHaloRowsFromOwners(geometry); } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + Jacobian.TransposeInPlace(); + } switch (KindPrecond) { case ILU: @@ -1700,6 +1704,8 @@ unsigned long CSysSolve::Solve_b(CSysMatrix& Jacobian, c /*--- If there was no call to solve first the preconditioner needs to be built here. ---*/ if (directCall) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { Jacobian.CopyHaloRowsFromOwners(geometry); } + END_SU2_OMP_SAFE_GLOBAL_ACCESS Jacobian.TransposeInPlace(); normal_prec->Build(); } diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 250f96e66c12..421d17cd9768 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -347,7 +347,7 @@ def main(): da_sp_pinArray_cht_2d_mf.cfg_dir = "incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d" da_sp_pinArray_cht_2d_mf.cfg_file = "configMaster.cfg" da_sp_pinArray_cht_2d_mf.test_iter = 100 - da_sp_pinArray_cht_2d_mf.test_vals = [-4.544489, -1.200740, -1.451677, -0.809551, -5.998322, -15.723654, -49.833480] + da_sp_pinArray_cht_2d_mf.test_vals = [-4.545587, -1.201859, -1.452976, -0.809884, -5.995415, -15.724051, 2.620874] da_sp_pinArray_cht_2d_mf.multizone = True test_list.append(da_sp_pinArray_cht_2d_mf) From 50639a4bfc7a7a64906683fd1c336c0001875751 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 18 Oct 2026 20:22:45 +0200 Subject: [PATCH 6/8] Move the halo row copy into TransposeInPlace and tidy comments TransposeInPlace now takes the geometry and copies the owner rows into the halo rows itself, so the two callers cannot forget it. Comments shortened and the new code run through clang-format. Co-Authored-By: Claude Sonnet 5.5 --- Common/include/linear_algebra/CSysMatrix.hpp | 15 ++++++++------- Common/src/linear_algebra/CSysMatrix.cpp | 13 +++++++------ Common/src/linear_algebra/CSysSolve.cpp | 10 ++-------- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 11 +++++------ 4 files changed, 22 insertions(+), 27 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 1adc0187b4ee..c95aaa10a70a 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -465,6 +465,12 @@ class CSysMatrix { return SU2_TYPE::GetValue(val); } + /*! + * \brief Overwrites the halo rows with the rows of the owning ranks, which hold the boundary condition modifications. + * \param[in] geometry - Geometry with the point-to-point communication pattern. + */ + void CopyHaloRowsFromOwners(const CGeometry* geometry); + /*! * \brief Calculates the matrix-vector product: product = matrix*vector * \param[in] matrix @@ -1293,16 +1299,11 @@ class CSysMatrix { */ void SetDiagonalAsColumnSum(); - /*! - * \brief Overwrites the halo rows with the rows of the owning ranks, which hold the boundary condition modifications. - * \param[in] geometry - Geometry with the point-to-point communication pattern. - */ - void CopyHaloRowsFromOwners(const CGeometry* geometry); - /*! * \brief Transposes the matrix, any preconditioner that was computed may be invalid. + * \param[in] geometry - Needed to give the halo rows the boundary condition modifications of their owners. */ - void TransposeInPlace(); + void TransposeInPlace(const CGeometry* geometry); /*! * \brief Add a scaled sparse matrix to "this" (axpy-type operation, A = A+alpha*B). diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 8a25e085cb41..b64c4993d717 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -1716,8 +1716,7 @@ void CSysMatrix::CopyHaloRowsFromOwners(const CGeometry* geometry) { const auto* nodes = geometry->nodes; const auto blkSize = nVar * nEqn; - /*--- Each point sent to a rank is followed by its blocks in the columns the sender does not own, - * as the block count, then the global column index and the block for each. ---*/ + /*--- Per sent point: block count, then global column and block for each column not owned by the sender. ---*/ vector> sendTo(size); for (int iMsg = 0; iMsg < geometry->nP2PSend; iMsg++) { auto& buf = sendTo[geometry->Neighbors_P2PSend[iMsg]]; @@ -1749,8 +1748,7 @@ void CSysMatrix::CopyHaloRowsFromOwners(const CGeometry* geometry) { MPIWrapper::Alltoallv(sendBuf.data(), nSend.data(), sendDisp.data(), MPI_DOUBLE, recvBuf.data(), nRecv.data(), recvDisp.data(), MPI_DOUBLE, SU2_MPI::GetComm()); - /*--- The points of a message arrive in the order of the receive list from that rank. Only the - * blocks in owned columns are copied, these become owned rows when the matrix is transposed. ---*/ + /*--- Points arrive in receive-list order, only blocks in columns owned here are kept. ---*/ for (int iMsg = 0; iMsg < geometry->nP2PRecv; iMsg++) { const auto* data = recvBuf.data() + recvDisp[geometry->Neighbors_P2PRecv[iMsg]]; for (auto k = geometry->nPoint_P2PRecv[iMsg]; k < geometry->nPoint_P2PRecv[iMsg + 1]; k++) { @@ -1761,7 +1759,8 @@ void CSysMatrix::CopyHaloRowsFromOwners(const CGeometry* geometry) { for (const auto jPoint : nodes->GetPoints(iPoint)) { if (!nodes->GetDomain(jPoint) || nodes->GetGlobalIndex(jPoint) != jGlobal) continue; auto* blk = GetBlock(iPoint, jPoint); - if (blk) for (auto i = 0ul; i < blkSize; i++) blk[i] = data[i]; + if (blk) + for (auto i = 0ul; i < blkSize; i++) blk[i] = data[i]; break; } } @@ -1770,8 +1769,10 @@ void CSysMatrix::CopyHaloRowsFromOwners(const CGeometry* geometry) { } template -void CSysMatrix::TransposeInPlace() { +void CSysMatrix::TransposeInPlace(const CGeometry* geometry) { SU2_ZONE_SCOPED + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { CopyHaloRowsFromOwners(geometry); } + END_SU2_OMP_SAFE_GLOBAL_ACCESS assert(nVar == nEqn && "Cannot transpose with nVar != nEqn."); auto swapAndTransp = [](unsigned long n, ScalarType* a, ScalarType* b) { diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 2ed0205f9a0b..a75758984ea8 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -1590,11 +1590,7 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con /*--- Build preconditioner for the transposed Jacobian ---*/ - if (RequiresTranspose) { - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { Jacobian.CopyHaloRowsFromOwners(geometry); } - END_SU2_OMP_SAFE_GLOBAL_ACCESS - Jacobian.TransposeInPlace(); - } + if (RequiresTranspose) Jacobian.TransposeInPlace(geometry); switch (KindPrecond) { case ILU: @@ -1704,9 +1700,7 @@ unsigned long CSysSolve::Solve_b(CSysMatrix& Jacobian, c /*--- If there was no call to solve first the preconditioner needs to be built here. ---*/ if (directCall) { - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { Jacobian.CopyHaloRowsFromOwners(geometry); } - END_SU2_OMP_SAFE_GLOBAL_ACCESS - Jacobian.TransposeInPlace(); + Jacobian.TransposeInPlace(geometry); normal_prec->Build(); } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index cc832823d7b8..9add40453933 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -3665,7 +3665,7 @@ unsigned long CIncEulerSolver::RegisterSolutionExtra(bool input, const CConfig* void CIncEulerSolver::SetAdjoint_SolutionExtra(const su2activevector& adj_sol, const CConfig* config) { SU2_ZONE_SCOPED if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { - /*--- Every rank holds a copy of the pressure drop, the reverse Allreduce of the massflow sums their seeds. ---*/ + /*--- The reverse Allreduce of the massflow sums the seeds of the copies on all ranks. ---*/ SU2_TYPE::SetDerivative(SPvalsUpdated.Streamwise_Periodic_PressureDrop, SU2_TYPE::GetValue(adj_sol[0]) / size); } } @@ -3675,11 +3675,10 @@ void CIncEulerSolver::ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, con if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { /*--- Sum the adjoints of the copies of the pressure drop on all ranks. ---*/ SU2_OMP_SAFE_GLOBAL_ACCESS( - const passivedouble local = SU2_TYPE::GetDerivative(SPvals.Streamwise_Periodic_PressureDrop); - passivedouble global = 0.0; - SelectMPIWrapper::W::Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - adj_sol[0] = global; - ) + const passivedouble local = SU2_TYPE::GetDerivative(SPvals.Streamwise_Periodic_PressureDrop); + passivedouble global = 0.0; + SelectMPIWrapper::W::Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + adj_sol[0] = global;) } } From 803684ffd3fd56ce52947f8b75139acaff2e853f Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 23 Oct 2026 15:02:13 +0200 Subject: [PATCH 7/8] Address review: integer halo indices and load-bounded line moves Send global column indices in an unsigned long buffer, apart from the scalar blocks, so single-precision builds keep indices above 2^24 exact. Reject whole-line moves in KeepLinesTogether when the resulting vertex weight imbalance exceeds max(1.5, the ParMETIS partition's imbalance). Add unit tests for the halo block buffer and the imbalance measure. Co-Authored-By: Claude Sonnet 5.5 --- Common/include/geometry/CPhysicalGeometry.hpp | 3 +- .../include/toolboxes/halo_block_buffer.hpp | 90 +++++++++++++++++++ Common/include/toolboxes/partition_load.hpp | 42 +++++++++ Common/src/geometry/CPhysicalGeometry.cpp | 24 +++-- Common/src/linear_algebra/CSysMatrix.cpp | 60 ++++++++----- .../Common/toolboxes/partition_tests.cpp | 62 +++++++++++++ UnitTests/meson.build | 3 +- 7 files changed, 251 insertions(+), 33 deletions(-) create mode 100644 Common/include/toolboxes/halo_block_buffer.hpp create mode 100644 Common/include/toolboxes/partition_load.hpp create mode 100644 UnitTests/Common/toolboxes/partition_tests.cpp diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp index e69b1a04a9b0..2994fe339f26 100644 --- a/Common/include/geometry/CPhysicalGeometry.hpp +++ b/Common/include/geometry/CPhysicalGeometry.hpp @@ -365,11 +365,12 @@ class CPhysicalGeometry final : public CGeometry { /*! * \brief Moves each line to the partition that holds most of its points. * \param[in] root - Line roots from LinkAnisotropicLines. + * \param[in] vwgt - Vertex weights ParMETIS balances; no move is made if it would unbalance them. * \param[in,out] part - Partition of every local point. * \param[out] longestLine - Number of points of the longest line. * \return Number of points that changed partition, summed over all ranks. */ - unsigned long KeepLinesTogether(const vector& root, vector& part, + unsigned long KeepLinesTogether(const vector& root, const vector& vwgt, vector& part, unsigned long& longestLine) const; #endif diff --git a/Common/include/toolboxes/halo_block_buffer.hpp b/Common/include/toolboxes/halo_block_buffer.hpp new file mode 100644 index 000000000000..622b46f002de --- /dev/null +++ b/Common/include/toolboxes/halo_block_buffer.hpp @@ -0,0 +1,90 @@ +/*! + * \file halo_block_buffer.hpp + * \brief Message buffers for matrix blocks, with integer indices kept apart from scalar values. + * \author SU2 Contributors + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include + +/*! + * \class CHaloBlockBuffer + * \brief Per point: a block count, then a global column index and a block per entry, in separate arrays. + */ +template +struct CHaloBlockBuffer { + std::vector index; + std::vector value; + + /*! + * \brief Starts a new point. + * \return Position of the point's block count. + */ + size_t BeginPoint() { + index.push_back(0); + return index.size() - 1; + } + + /*! + * \brief Appends a block of n scalars to the point started at countPos. + */ + void AddBlock(size_t countPos, unsigned long column, const ScalarType* block, size_t n) { + index.push_back(column); + value.insert(value.end(), block, block + n); + index[countPos] += 1; + } +}; + +/*! + * \class CHaloBlockReader + * \brief Walks the arrays written by CHaloBlockBuffer. + */ +template +class CHaloBlockReader { + const unsigned long* idx; + const ScalarType* val; + size_t blkSize; + + public: + CHaloBlockReader(const unsigned long* index, const ScalarType* value, size_t blockSize) + : idx(index), val(value), blkSize(blockSize) {} + + /*! + * \return Number of blocks of the next point. + */ + unsigned long NextPoint() { return *idx++; } + + /*! + * \brief Reads the next block of the current point. + * \param[out] column - Its global column index. + * \return Pointer to its scalars. + */ + const ScalarType* NextBlock(unsigned long& column) { + column = *idx++; + const auto* blk = val; + val += blkSize; + return blk; + } +}; diff --git a/Common/include/toolboxes/partition_load.hpp b/Common/include/toolboxes/partition_load.hpp new file mode 100644 index 000000000000..5cdedd26759e --- /dev/null +++ b/Common/include/toolboxes/partition_load.hpp @@ -0,0 +1,42 @@ +/*! + * \file partition_load.hpp + * \brief Load balance measure for graph partitions. + * \author SU2 Contributors + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include +#include +#include + +/*! + * \brief Ratio of the largest partition load to the mean load. + * \param[in] load - Total vertex weight of each partition. + */ +inline double PartitionImbalance(const std::vector& load) { + const auto sum = std::accumulate(load.begin(), load.end(), 0ul); + if (sum == 0) return 1.0; + return static_cast(*std::max_element(load.begin(), load.end())) * load.size() / sum; +} diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 382e5a6e5748..bcb48e0d5ccf 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -29,6 +29,7 @@ #include "../../include/adt/CADTPointsOnlyClass.hpp" #include "../../include/toolboxes/printing_toolbox.hpp" #include "../../include/toolboxes/CLinearPartitioner.hpp" +#include "../../include/toolboxes/partition_load.hpp" #include "../../include/toolboxes/C1DInterpolation.hpp" #include "../../include/toolboxes/geometry_toolbox.hpp" #include "../../include/geometry/meshreader/CSU2ASCIIMeshReaderFEM.hpp" @@ -7449,8 +7450,8 @@ vector CPhysicalGeometry::LinkAnisotropicLines(const vector& root, vector& part, - unsigned long& longestLine) const { +unsigned long CPhysicalGeometry::KeepLinesTogether(const vector& root, const vector& vwgt, + vector& part, unsigned long& longestLine) const { CLinearPartitioner pointPartitioner(Global_nPointDomain, 0); const unsigned long firstIndex = pointPartitioner.GetFirstIndexOnRank(rank); @@ -7496,12 +7497,19 @@ unsigned long CPhysicalGeometry::KeepLinesTogether(const vector& decision[line.first - firstIndex] = best; } - /*--- Every point takes the partition of its line, unless that would leave a rank without points. ---*/ + /*--- Every point takes the partition of its line, unless that leaves the load more unbalanced than + * both the limit and the partition ParMETIS returned. ---*/ const auto chosen = FetchFromOwners(root, decision); - vector count(size, 0), countGlobal(size); - for (const auto iPart : chosen) count[iPart]++; - SU2_MPI::Allreduce(count.data(), countGlobal.data(), size, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); - if (*min_element(countGlobal.begin(), countGlobal.end()) == 0) return 0; + vector oldLoad(size, 0), newLoad(size, 0), oldLoadGlobal(size), newLoadGlobal(size); + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + oldLoad[part[iPoint]] += vwgt[iPoint]; + newLoad[chosen[iPoint]] += vwgt[iPoint]; + } + SU2_MPI::Allreduce(oldLoad.data(), oldLoadGlobal.data(), size, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(newLoad.data(), newLoadGlobal.data(), size, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + constexpr double MAX_LINE_IMBALANCE = 1.5; + if (*min_element(newLoadGlobal.begin(), newLoadGlobal.end()) == 0) return 0; + if (PartitionImbalance(newLoadGlobal) > max(MAX_LINE_IMBALANCE, PartitionImbalance(oldLoadGlobal))) return 0; unsigned long moved = 0, movedGlobal = 0; for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { @@ -7722,7 +7730,7 @@ void CPhysicalGeometry::SetColorGrid_Parallel(const CConfig* config) { if (anisotropic) { unsigned long longestLine = 0; - const auto moved = KeepLinesTogether(lineRoot, part, longestLine); + const auto moved = KeepLinesTogether(lineRoot, vwgt, part, longestLine); /*--- The imbalance is measured in the vertex weights ParMETIS balances. ---*/ vector load(size, 0), loadGlobal(size); for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) load[part[iPoint]] += vwgt[iPoint]; diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index b64c4993d717..0f163320d1a9 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -30,6 +30,7 @@ #include "../../include/geometry/CGeometry.hpp" #include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/toolboxes/allocation_toolbox.hpp" +#include "../../include/toolboxes/halo_block_buffer.hpp" #include #include @@ -1716,46 +1717,59 @@ void CSysMatrix::CopyHaloRowsFromOwners(const CGeometry* geometry) { const auto* nodes = geometry->nodes; const auto blkSize = nVar * nEqn; - /*--- Per sent point: block count, then global column and block for each column not owned by the sender. ---*/ - vector> sendTo(size); + /*--- Per sent point: block count, then global column and block for each column not owned by the sender. + * Indices and coefficients travel in separate buffers so that indices keep full integer precision. ---*/ + vector> sendTo(size); + vector tmp(blkSize); for (int iMsg = 0; iMsg < geometry->nP2PSend; iMsg++) { auto& buf = sendTo[geometry->Neighbors_P2PSend[iMsg]]; for (auto k = geometry->nPoint_P2PSend[iMsg]; k < geometry->nPoint_P2PSend[iMsg + 1]; k++) { const auto iPoint = geometry->Local_Point_P2PSend[k]; - const auto countPos = buf.size(); - buf.push_back(0); + const auto countPos = buf.BeginPoint(); for (const auto jPoint : nodes->GetPoints(iPoint)) { if (nodes->GetDomain(jPoint)) continue; const auto* blk = GetBlock(iPoint, jPoint); if (!blk) continue; - buf.push_back(nodes->GetGlobalIndex(jPoint)); - for (auto i = 0ul; i < blkSize; i++) buf.push_back(SU2_TYPE::GetValue(blk[i])); - buf[countPos] += 1; + for (auto i = 0ul; i < blkSize; i++) tmp[i] = SU2_TYPE::GetValue(blk[i]); + buf.AddBlock(countPos, nodes->GetGlobalIndex(jPoint), tmp.data(), blkSize); } } } - vector nSend(size), nRecv(size), sendDisp(size + 1, 0), recvDisp(size + 1, 0); - for (int iRank = 0; iRank < size; iRank++) nSend[iRank] = sendTo[iRank].size(); - MPIWrapper::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); - for (int iRank = 0; iRank < size; iRank++) { - sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; - recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; - } - vector sendBuf(sendDisp[size]), recvBuf(recvDisp[size]); - for (int iRank = 0; iRank < size; iRank++) - copy(sendTo[iRank].begin(), sendTo[iRank].end(), sendBuf.begin() + sendDisp[iRank]); - MPIWrapper::Alltoallv(sendBuf.data(), nSend.data(), sendDisp.data(), MPI_DOUBLE, recvBuf.data(), nRecv.data(), - recvDisp.data(), MPI_DOUBLE, SU2_MPI::GetComm()); + /*--- Exchanges one array per rank, returns the received data and the per-rank displacements. ---*/ + auto exchange = [&](auto getArray, auto dummy, auto mpiType, auto& recvBuf, vector& recvDisp) { + using T = decltype(dummy); + vector nSend(size), nRecv(size), sendDisp(size + 1, 0); + recvDisp.assign(size + 1, 0); + for (int iRank = 0; iRank < size; iRank++) nSend[iRank] = getArray(sendTo[iRank]).size(); + MPIWrapper::Alltoall(nSend.data(), 1, MPI_INT, nRecv.data(), 1, MPI_INT, SU2_MPI::GetComm()); + for (int iRank = 0; iRank < size; iRank++) { + sendDisp[iRank + 1] = sendDisp[iRank] + nSend[iRank]; + recvDisp[iRank + 1] = recvDisp[iRank] + nRecv[iRank]; + } + vector sendBuf(sendDisp[size]); + recvBuf.resize(recvDisp[size]); + for (int iRank = 0; iRank < size; iRank++) + copy(getArray(sendTo[iRank]).begin(), getArray(sendTo[iRank]).end(), sendBuf.begin() + sendDisp[iRank]); + MPIWrapper::Alltoallv(sendBuf.data(), nSend.data(), sendDisp.data(), mpiType, recvBuf.data(), nRecv.data(), + recvDisp.data(), mpiType, SU2_MPI::GetComm()); + }; + vector recvIdx; + vector recvVal; + vector idxDisp, valDisp; + exchange([](auto& b) -> auto& { return b.index; }, 0ul, MPI_UNSIGNED_LONG, recvIdx, idxDisp); + exchange([](auto& b) -> auto& { return b.value; }, passivedouble(0), MPI_DOUBLE, recvVal, valDisp); /*--- Points arrive in receive-list order, only blocks in columns owned here are kept. ---*/ for (int iMsg = 0; iMsg < geometry->nP2PRecv; iMsg++) { - const auto* data = recvBuf.data() + recvDisp[geometry->Neighbors_P2PRecv[iMsg]]; + const auto iRank = geometry->Neighbors_P2PRecv[iMsg]; + CHaloBlockReader reader(recvIdx.data() + idxDisp[iRank], recvVal.data() + valDisp[iRank], blkSize); for (auto k = geometry->nPoint_P2PRecv[iMsg]; k < geometry->nPoint_P2PRecv[iMsg + 1]; k++) { const auto iPoint = geometry->Local_Point_P2PRecv[k]; - const auto nBlk = static_cast(*data++); - for (auto iBlk = 0ul; iBlk < nBlk; iBlk++, data += blkSize) { - const auto jGlobal = static_cast(*data++); + const auto nBlk = reader.NextPoint(); + for (auto iBlk = 0ul; iBlk < nBlk; iBlk++) { + unsigned long jGlobal; + const auto* data = reader.NextBlock(jGlobal); for (const auto jPoint : nodes->GetPoints(iPoint)) { if (!nodes->GetDomain(jPoint) || nodes->GetGlobalIndex(jPoint) != jGlobal) continue; auto* blk = GetBlock(iPoint, jPoint); diff --git a/UnitTests/Common/toolboxes/partition_tests.cpp b/UnitTests/Common/toolboxes/partition_tests.cpp new file mode 100644 index 000000000000..b18244446719 --- /dev/null +++ b/UnitTests/Common/toolboxes/partition_tests.cpp @@ -0,0 +1,62 @@ +/*! + * \file partition_tests.cpp + * \brief Unit tests for the halo block message buffers and the partition imbalance measure. + * \author SU2 Contributors + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include "../../../Common/include/toolboxes/halo_block_buffer.hpp" +#include "../../../Common/include/toolboxes/partition_load.hpp" + +TEST_CASE("Halo block buffer keeps neighboring large indices apart", "[Toolboxes]") { + /*--- 2^24 and 2^24 + 1 are the same float. ---*/ + constexpr unsigned long big = 16777216; + const float blockA[2] = {1.f, 2.f}, blockB[2] = {3.f, 4.f}; + + CHaloBlockBuffer buf; + const auto pos = buf.BeginPoint(); + buf.AddBlock(pos, big, blockA, 2); + buf.AddBlock(pos, big + 1, blockB, 2); + const auto pos2 = buf.BeginPoint(); + buf.AddBlock(pos2, 7, blockB, 2); + + CHaloBlockReader reader(buf.index.data(), buf.value.data(), 2); + unsigned long col; + REQUIRE(reader.NextPoint() == 2); + const auto* a = reader.NextBlock(col); + CHECK(col == big); + CHECK(a[0] == 1.f); + const auto* b = reader.NextBlock(col); + CHECK(col == big + 1); + CHECK(b[0] == 3.f); + REQUIRE(reader.NextPoint() == 1); + reader.NextBlock(col); + CHECK(col == 7); +} + +TEST_CASE("Partition imbalance flags moving a dominant line", "[Toolboxes]") { + CHECK(PartitionImbalance({255, 255, 255, 255}) == Approx(1.0)); + CHECK(PartitionImbalance({930, 30, 30, 30}) == Approx(3.64706).epsilon(1e-4)); + CHECK(PartitionImbalance({0, 0}) == Approx(1.0)); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index c6216606338f..f2bd69a0be02 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -25,7 +25,8 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'SU2_CFD/windowing.cpp', 'Common/toolboxes/random_toolbox_tests.cpp', 'Common/linear_algebra/quantization_tests.cpp', - 'Common/linear_algebra/edge_residual_blocks_tests.cpp']) + 'Common/linear_algebra/edge_residual_blocks_tests.cpp', + 'Common/toolboxes/partition_tests.cpp']) # Reverse-mode (algorithmic differentiation) tests: su2_cfd_tests_ad = files(['Common/simple_ad_test.cpp', From 51cadb3dba8a0f0930920128607ae51a0f99bb03 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 23 Oct 2026 15:16:36 +0200 Subject: [PATCH 8/8] pre-commit --- Common/src/linear_algebra/CSysMatrix.cpp | 6 ++++-- TestCases/vandv.py | 18 +++++++++--------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 0f163320d1a9..eaf2b9a3776a 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -1757,8 +1757,10 @@ void CSysMatrix::CopyHaloRowsFromOwners(const CGeometry* geometry) { vector recvIdx; vector recvVal; vector idxDisp, valDisp; - exchange([](auto& b) -> auto& { return b.index; }, 0ul, MPI_UNSIGNED_LONG, recvIdx, idxDisp); - exchange([](auto& b) -> auto& { return b.value; }, passivedouble(0), MPI_DOUBLE, recvVal, valDisp); + exchange( + [](auto& b) -> auto& { return b.index; }, 0ul, MPI_UNSIGNED_LONG, recvIdx, idxDisp); + exchange( + [](auto& b) -> auto& { return b.value; }, passivedouble(0), MPI_DOUBLE, recvVal, valDisp); /*--- Points arrive in receive-list order, only blocks in columns owned here are kept. ---*/ for (int iMsg = 0; iMsg < geometry->nP2PRecv; iMsg++) { diff --git a/TestCases/vandv.py b/TestCases/vandv.py index 88602b2d2275..b2b2506f02ae 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -45,7 +45,7 @@ def main(): p30n30.cfg_dir = "vandv/rans/30p30n" p30n30.cfg_file = "config.cfg" p30n30.test_iter = 5 - p30n30.test_vals = [-11.267106, -11.168215, -11.182822, -10.949673, -14.233489, 0.052235, 2.830394, 1.318894, -1.210646, 1.000000, 12.763000] + p30n30.test_vals = [-11.251502, -11.158987, -11.170194, -10.927262, -14.233490, 0.052235, 2.830394, 1.318894, -1.219937, 1.000000, 12.763000] test_list.append(p30n30) # This is not part of the V&V cases yet, its tested in this script because it is a relatively long test (~1 min). @@ -54,7 +54,7 @@ def main(): p30n30_ad.cfg_dir = "vandv/rans/30p30n" p30n30_ad.cfg_file = "config_ad.cfg" p30n30_ad.test_iter = 9 - p30n30_ad.test_vals = [-7.283709, -6.072615, -5.995304, -7.197048, -4.568373, -1.167146, -2.316777, 1.1791e-01, 3.4123e+01] + p30n30_ad.test_vals = [-7.205798, -5.937759, -5.846810, -7.076967, -4.534548, -1.226612, -2.281655, 0.117910, 34.123000] p30n30_ad.command = TestCase.Command("mpirun -n 2", "SU2_CFD_AD") test_list.append(p30n30_ad) @@ -63,7 +63,7 @@ def main(): flatplate_sst1994m.cfg_dir = "vandv/rans/flatplate" flatplate_sst1994m.cfg_file = "turb_flatplate_sst.cfg" flatplate_sst1994m.test_iter = 5 - flatplate_sst1994m.test_vals = [-13.040943, -10.136971, -10.942764, -7.985181, -10.323857, -4.732487, 0.002801] + flatplate_sst1994m.test_vals = [-13.041850, -10.137145, -10.934937, -7.992806, -10.323854, -4.732544, 0.002801] flatplate_sst1994m.test_vals_aarch64 = [-13.021715, -9.534786, -10.401912, -7.501836, -9.750800, -4.850665, 0.002807] test_list.append(flatplate_sst1994m) @@ -72,7 +72,7 @@ def main(): bump_sst1994m.cfg_dir = "vandv/rans/bump_in_channel" bump_sst1994m.cfg_file = "turb_bump_sst.cfg" bump_sst1994m.test_iter = 5 - bump_sst1994m.test_vals = [-11.928292, -10.095796, -9.512953, -6.445652, -11.774088, -6.988752, 0.004931] + bump_sst1994m.test_vals = [-11.927778, -10.095211, -9.513000, -6.445108, -11.774670, -6.990602, 0.004931] bump_sst1994m.test_vals_aarch64 = [-13.042689, -10.812982, -10.604523, -7.655547, -10.816257, -5.308083, 0.004911] test_list.append(bump_sst1994m) @@ -81,7 +81,7 @@ def main(): swbli_sa.cfg_dir = "vandv/rans/swbli" swbli_sa.cfg_file = "config_sa.cfg" swbli_sa.test_iter = 5 - swbli_sa.test_vals = [-11.504424, -10.941741, -12.049925, -10.586263, -16.090385, 0.002242, -1.614365, 1.340100] + swbli_sa.test_vals = [-11.504424, -10.941741, -12.049926, -10.586263, -16.090438, 0.002242, -1.614397, 1.34010] swbli_sa.test_vals_aarch64 = [-11.504424, -10.941741, -12.049925, -10.586263, -16.090385, 0.002242, -1.614365, 1.340100] test_list.append(swbli_sa) @@ -91,7 +91,7 @@ def main(): swbli_sst.cfg_dir = "vandv/rans/swbli" swbli_sst.cfg_file = "config_sst.cfg" swbli_sst.test_iter = 5 - swbli_sst.test_vals = [-11.569218, -10.909086, -11.607984, -10.431163, -11.407588, -2.637660, 0.001816, -1.305818, -3.514590, 13.399000] + swbli_sst.test_vals = [-11.567956, -10.908257, -11.607992, -10.430471, -11.407660, -2.637661, 0.001816, -1.354761, -3.514714, 13.399000] test_list.append(swbli_sst) # DSMA661 - SA @@ -99,7 +99,7 @@ def main(): dsma661_sa.cfg_dir = "vandv/rans/dsma661" dsma661_sa.cfg_file = "dsma661_sa_config.cfg" dsma661_sa.test_iter = 5 - dsma661_sa.test_vals = [-11.255379, -8.243275, -9.037314, -5.927070, -10.737675, 0.155687, 0.024232] + dsma661_sa.test_vals = [-11.247167, -8.241024, -9.035101, -5.907729, -10.737682, 0.155687, 0.024232] dsma661_sa.test_vals_aarch64 = [-11.293183, -8.241775, -9.083761, -6.011398, -10.737680, 0.155687, 0.024232] test_list.append(dsma661_sa) @@ -108,7 +108,7 @@ def main(): dsma661_sst.cfg_dir = "vandv/rans/dsma661" dsma661_sst.cfg_file = "dsma661_sst_config.cfg" dsma661_sst.test_iter = 5 - dsma661_sst.test_vals = [-11.025153, -8.156995, -9.057021, -5.947228, -10.650874, -7.884423, 0.155882, 0.023344] + dsma661_sst.test_vals = [-11.023380, -8.156126, -9.055628, -5.935182, -10.651455, -7.911250, 0.155882, 0.023344] dsma661_sst.test_vals_aarch64 = [-10.977195, -8.403731, -8.747068, -5.808899, -10.522786, -7.369851, 0.155875, 0.023353] test_list.append(dsma661_sst) @@ -121,7 +121,7 @@ def main(): sandiajet_sst.cfg_dir = "vandv/species_transport/sandia_jet" sandiajet_sst.cfg_file = "validation.cfg" sandiajet_sst.test_iter = 5 - sandiajet_sst.test_vals = [-4.930056, -3.072757, -3.531917, -3.632319, 1.186060, -6.175518, 5.000000, -1.413277, 5.000000, -4.735257, 5.000000, -1.924560, 0.000256, 0.000000, 0.000000, 0.000256, 4126.200000, 4008.100000, 62.399000, 55.713000] + sandiajet_sst.test_vals = [-4.922130, -3.073873, -3.517555, -3.626843, 1.094593, -6.182863, 5.000000, -1.310699, 4.000000, -5.581570, 5.000000, -3.793476, 0.000256, 0.000000, 0.000000, 0.000256, 4127.300000, 4009.300000, 62.402000, 55.615000] sandiajet_sst.test_vals_aarch64 = [-8.271400, -6.305945, -6.123143, -5.634565, -1.687556, -7.869775, 5.000000, -1.527449, 5.000000, -4.053788, 5.000000, -2.303502, 0.000257, 0.000000, 0.000000, 0.000257, 4020.500000, 3919.900000, 49.151000, 51.436000] test_list.append(sandiajet_sst)