Partition anisotropic lines onto one rank and fix adjoint mpi ranks

.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: |

Common/include/CConfig.hpp

@@ -1117,6 +1117,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. */
@@ -10334,6 +10335,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.

Common/include/geometry/CPhysicalGeometry.hpp

@@ -337,6 +337,43 @@ 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<idx_t> 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<unsigned long> FetchFromOwners(const vector<unsigned long>& wanted,
+                                        const vector<unsigned long>& 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<unsigned long> LinkAnisotropicLines(const vector<idx_t>& adjwgt, vector<char>& lineEdge) const;
+
+  /*!
+   * \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<unsigned long>& root, const vector<idx_t>& vwgt, vector<idx_t>& 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.

Common/include/linear_algebra/CSysMatrix.hpp

@@ -468,6 +468,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
@@ -1298,8 +1304,9 @@ class CSysMatrix {
 
   /*!
    * \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).

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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <vector>
+
+/*!
+ * \class CHaloBlockBuffer
+ * \brief Per point: a block count, then a global column index and a block per entry, in separate arrays.
+ */
+template <class ScalarType>
+struct CHaloBlockBuffer {
+  std::vector<unsigned long> index;
+  std::vector<ScalarType> 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 ScalarType>
+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;
+  }
+};

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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <numeric>
+#include <vector>
+
+/*!
+ * \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<unsigned long>& load) {
+  const auto sum = std::accumulate(load.begin(), load.end(), 0ul);
+  if (sum == 0) return 1.0;
+  return static_cast<double>(*std::max_element(load.begin(), load.end())) * load.size() / sum;
+}

Common/src/CConfig.cpp

@@ -3115,6 +3115,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*/
 
@@ -4889,6 +4892,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;