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/include/CConfig.hpp b/Common/include/CConfig.hpp
index dd393924bb36..71e39296c2de 100644
--- a/Common/include/CConfig.hpp
+++ b/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.
diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp
index 832dce0e1249..2e13b471173f 100644
--- a/Common/include/geometry/CPhysicalGeometry.hpp
+++ b/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.
diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp
index ed33589b31d2..67815cfd7fb0 100644
--- a/Common/include/linear_algebra/CSysMatrix.hpp
+++ b/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).
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 <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;
+  }
+};
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 <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;
+}
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index f2b171c9571f..deb3860549d8 100644
--- a/Common/src/CConfig.cpp
+++ b/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;
diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp
index fa9d6b42e3fb..9114c11fd99b 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"
@@ -2426,9 +2427,9 @@ void CPhysicalGeometry::LoadSurfaceElements(CConfig* config, CGeometry* geometry
 
   unsigned long Local_Nodes[N_POINTS_HEXAHEDRON];
 
-  vector<vector<unsigned long> > Line_List;
-  vector<vector<unsigned long> > BoundTria_List;
-  vector<vector<unsigned long> > BoundQuad_List;
+  vector<vector<unsigned long>> Line_List;
+  vector<vector<unsigned long>> BoundTria_List;
+  vector<vector<unsigned long>> BoundQuad_List;
 
   vector<unsigned long> Marker_Local;
 
@@ -2967,14 +2968,14 @@ void CPhysicalGeometry::SetSendReceive(const CConfig* config) {
   unsigned short iNode, jNode;
   vector<unsigned long>::iterator it;
 
-  vector<vector<unsigned long> >
+  vector<vector<unsigned long>>
       SendTransfLocal; /*!< \brief Vector to store the type of transformation for this send point. */
-  vector<vector<unsigned long> >
+  vector<vector<unsigned long>>
       ReceivedTransfLocal; /*!< \brief Vector to store the type of transformation for this received point. */
-  vector<vector<unsigned long> > SendDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the point
-                                                     index of the node that must me sended. */
-  vector<vector<unsigned long> > ReceivedDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the
-                                                         point index of the node that must me sended. */
+  vector<vector<unsigned long>> SendDomainLocal;     /*!< \brief SendDomain[from domain][to domain] and return the point
+                                                        index of the node that must me sended. */
+  vector<vector<unsigned long>> ReceivedDomainLocal; /*!< \brief SendDomain[from domain][to domain] and return the
+                                                        point index of the node that must me sended. */
 
   unordered_map<unsigned long, unsigned long>::const_iterator MI;
 
@@ -4426,7 +4427,7 @@ void CPhysicalGeometry::SetPositive_ZArea(CConfig* config) {
 }
 
 void CPhysicalGeometry::SetPoint_Connectivity() {
-  vector<vector<unsigned long> > points(nPoint);
+  vector<vector<unsigned long>> points(nPoint);
 
   SU2_OMP_PARALLEL {
     unsigned short Node_Neighbor, iNode, iNeighbor;
@@ -4434,7 +4435,7 @@ void CPhysicalGeometry::SetPoint_Connectivity() {
 
     /*--- Loop over all the elements ---*/
     BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
-      vector<vector<long> > elems(nPoint);
+      vector<vector<long>> elems(nPoint);
 
       for (iElem = 0; iElem < nElem; iElem++) {
         /*--- Loop over all the nodes of an element ---*/
@@ -7361,6 +7362,285 @@ void CPhysicalGeometry::SetBoundTecPlot(char mesh_filename[MAX_STRING_SIZE], boo
   Tecplot_File.close();
 }
 
+#if defined(HAVE_MPI) && defined(HAVE_PARMETIS)
+vector<unsigned long> CPhysicalGeometry::FetchFromOwners(const vector<unsigned long>& wanted,
+                                                         const vector<unsigned long>& localValue) const {
+  CLinearPartitioner pointPartitioner(Global_nPointDomain, 0);
+  const unsigned long firstIndex = pointPartitioner.GetFirstIndexOnRank(rank);
+
+  vector<int> 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<unsigned long> sendIndex(sendDisp[size]), slotOf(wanted.size()), recvIndex(recvDisp[size]);
+  vector<int> 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<unsigned long> 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<unsigned long> CPhysicalGeometry::LinkAnisotropicLines(const vector<idx_t>& adjwgt,
+                                                              vector<char>& 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<unsigned long> neighbor(adjacency.begin(), adjacency.end());
+
+  vector<unsigned long> 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<unsigned long>(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<unsigned long>(adjwgt[k]);
+      if (w < static_cast<unsigned long>(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<unsigned long> 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<unsigned long>& root, const vector<idx_t>& vwgt,
+                                                   vector<idx_t>& 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<int> 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<unsigned long> sendBuf(sendDisp[size]), recvBuf(recvDisp[size]);
+  vector<int> 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<unsigned long>(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<unsigned long, map<unsigned long, unsigned long>> 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<unsigned long> decision(nPoint);
+  for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) decision[iPoint] = static_cast<unsigned long>(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 leaves the load more unbalanced than
+   *    both the limit and the partition ParMETIS returned. ---*/
+  const auto chosen = FetchFromOwners(root, decision);
+  vector<unsigned long> 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++) {
+    if (static_cast<unsigned long>(part[iPoint]) != chosen[iPoint]) {
+      part[iPoint] = static_cast<idx_t>(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<idx_t> 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<vector<unsigned long>> 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<int> 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<int>(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<unsigned long> 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());
+
+  /*--- MPI_DOUBLE is an active type in AD builds, so the buffers hold su2double. ---*/
+  vector<su2double> 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] = 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<unsigned long, unsigned long> 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 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. ---*/
+
+  vector<idx_t> 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<idx_t>(elem[iElem]->GetNode(b));
+        const auto it = lower_bound(begin, end, jGlobal);
+        if ((it == end) || (*it != jGlobal)) continue;
+
+        const auto ratio = static_cast<idx_t>(sqrt(maxDist2 / d2));
+        auto& w = weight[it - adjacency.begin()];
+        w = max(w, min<idx_t>(maxWeight, max<idx_t>(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 +7692,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<idx_t> adjwgt;
+  vector<unsigned long> lineRoot;
+  if (anisotropic) {
+    constexpr long maxWeight = (sizeof(idx_t) >= 8) ? 1000000 : 1000;
+    adjwgt = ComputeAnisotropyEdgeWeights(maxWeight);
+    vector<char> 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 +7717,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, vwgt, part, longestLine);
+    /*--- The imbalance is measured in the vertex weights ParMETIS balances. ---*/
+    vector<unsigned long> 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/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp
index 1f46844c4918..56232a86714e 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 <cmath>
 #include <limits>
@@ -1707,8 +1708,87 @@ void CSysMatrix<ScalarType>::SetDiagonalAsColumnSum() {
 }
 
 template <class ScalarType>
-void CSysMatrix<ScalarType>::TransposeInPlace() {
+void CSysMatrix<ScalarType>::CopyHaloRowsFromOwners(const CGeometry* geometry) {
   SU2_ZONE_SCOPED
+  const int size = SU2_MPI::GetSize();
+  if (size == 1) return;
+
+  using MPIWrapper = SelectMPIWrapper<passivedouble>::W;
+  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.
+   *    Indices and coefficients travel in separate buffers so that indices keep full integer precision. ---*/
+  vector<CHaloBlockBuffer<passivedouble>> sendTo(size);
+  vector<passivedouble> 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.BeginPoint();
+      for (const auto jPoint : nodes->GetPoints(iPoint)) {
+        if (nodes->GetDomain(jPoint)) continue;
+        const auto* blk = GetBlock(iPoint, jPoint);
+        if (!blk) continue;
+        for (auto i = 0ul; i < blkSize; i++) tmp[i] = SU2_TYPE::GetValue(blk[i]);
+        buf.AddBlock(countPos, nodes->GetGlobalIndex(jPoint), tmp.data(), blkSize);
+      }
+    }
+  }
+
+  /*--- 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<int>& recvDisp) {
+    using T = decltype(dummy);
+    vector<int> 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<T> 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<unsigned long> recvIdx;
+  vector<passivedouble> recvVal;
+  vector<int> 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 iRank = geometry->Neighbors_P2PRecv[iMsg];
+    CHaloBlockReader<passivedouble> 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 = 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);
+          if (blk)
+            for (auto i = 0ul; i < blkSize; i++) blk[i] = data[i];
+          break;
+        }
+      }
+    }
+  }
+}
+
+template <class ScalarType>
+void CSysMatrix<ScalarType>::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 a2fc481c9784..5a37739bcf2f 100644
--- a/Common/src/linear_algebra/CSysSolve.cpp
+++ b/Common/src/linear_algebra/CSysSolve.cpp
@@ -1590,7 +1590,7 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
 
       /*--- Build preconditioner for the transposed Jacobian ---*/
 
-      if (RequiresTranspose) Jacobian.TransposeInPlace();
+      if (RequiresTranspose) Jacobian.TransposeInPlace(geometry);
 
       switch (KindPrecond) {
         case ILU:
@@ -1700,7 +1700,7 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
 
   /*--- If there was no call to solve first the preconditioner needs to be built here. ---*/
   if (directCall) {
-    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 08af0f78e72b..f11919329c30 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -3714,14 +3714,20 @@ 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]));
+    /*--- 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);
   }
 }
 
 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<passivedouble>::W::Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+        adj_sol[0] = global;)
   }
 }
 
diff --git a/SU2_PY/FSI_tools/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py
index b4df745e494e..a0b31bb193b2 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 f0def53fe935..93ad4fbda6a4 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 0c16c3a91de0..421ab049ff76 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 aa49d8a2788a..97ac22e565ac 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 55264dc5ff1c..10775e74b487 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 5fc61024aaab..3eb840aecb5d 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 5bc6247ddcfe..4e8216a68cf4 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 a7c7ea18d0f9..4071e2e9a3af 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 73ff43340473..07b9fd1ca6b3 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 94351d3d406e..00b7c2e5cacd 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 434ce0da4597..cd2f2e0fbe42 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 66a604af985a..4a978067524b 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 79282887bccd..ca8ddaf2f04a 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/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index c15cde2b2099..d2276a45e20d 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)
 
diff --git a/TestCases/vandv.py b/TestCases/vandv.py
index a2230114c32a..f67f2098c322 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.040942, -10.136971, -10.942969, -7.985184, -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.928283, -10.095826, -9.512968, -6.445644, -11.774076, -6.988223, 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.632034, 1.186060, -6.175515, 5.000000, -1.413277, 5.000000, -4.736007, 5.000000, -1.924583, 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)
 
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 <http://www.gnu.org/licenses/>.
+ */
+
+#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<float> 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<float> 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 2c9acd98b03d..fc932afc05e9 100644
--- a/UnitTests/meson.build
+++ b/UnitTests/meson.build
@@ -26,7 +26,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',
diff --git a/config_template.cfg b/config_template.cfg
index 66d7b191879b..073984da4c5d 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -2559,6 +2559,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)
diff --git a/meson.py b/meson.py
index b9ebcf644131..c879db5b5dd9 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")