diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp
index ca8a3a26aabd..a12ec045164e 100644
--- a/Common/include/geometry/CGeometry.hpp
+++ b/Common/include/geometry/CGeometry.hpp
@@ -319,6 +319,9 @@ class CGeometry {
   PassiveRequest* reqP_P2PSend{nullptr}; /*!< \brief Data structure for point-to-point send requests. */
   PassiveRequest* reqP_P2PRecv{nullptr}; /*!< \brief Data structure for point-to-point recv requests. */
 
+  bool periodicPairIntersections{false};    /*!< Whether distinct periodic pairs intersect in this zone. */
+  vector<passivedouble> periodicEdgeWeight; /*!< Fraction of each edge stencil owned by this periodic copy. */
+
   /*--- Data structures for periodic communications. ---*/
 
   int maxCountPerPeriodicPoint{0}; /*!< \brief Maximum number of pieces of data sent per vertex in periodic comms. */
@@ -346,8 +349,6 @@ class CGeometry {
                                               a particular vertex to be received in periodic comms. */
   su2double* bufD_PeriodicRecv{nullptr};   /*!< \brief Data structure for su2double periodic receive. */
   su2double* bufD_PeriodicSend{nullptr};   /*!< \brief Data structure for su2double periodic send. */
-  unsigned short* bufS_PeriodicRecv{nullptr};  /*!< \brief Data structure for unsigned long periodic receive. */
-  unsigned short* bufS_PeriodicSend{nullptr};  /*!< \brief Data structure for unsigned long periodic send. */
   SU2_MPI::Request* req_PeriodicSend{nullptr}; /*!< \brief Data structure for periodic send requests. */
   SU2_MPI::Request* req_PeriodicRecv{nullptr}; /*!< \brief Data structure for periodic recv requests. */
 
@@ -519,6 +520,18 @@ class CGeometry {
    */
   void AllocatePeriodicComms(unsigned short val_countPerPeriodicPoint);
 
+  /*! \brief Rebuild after changing connectivity, periodic membership or transforms; not needed for coordinate motion.
+   */
+  void PreprocessPeriodicGeometry(const CConfig& config);
+
+  /*! \brief Zone-wide decision, identical on all MPI ranks, for intersecting periodic stencils. */
+  inline bool HasPeriodicPairIntersections() const { return periodicPairIntersections; }
+
+  /*! \brief Fraction of an edge stencil owned by this periodic copy. */
+  inline passivedouble GetPeriodicEdgeWeight(unsigned long iEdge) const {
+    return periodicEdgeWeight.empty() ? 1.0 : periodicEdgeWeight[iEdge];
+  }
+
   /*!
    * \brief Routine to launch non-blocking recvs only for all periodic communication with neighboring partitions.
    * \note This routine is called by any class that has loaded data into the generic communication buffers.
diff --git a/Common/include/toolboxes/geometry_toolbox.hpp b/Common/include/toolboxes/geometry_toolbox.hpp
index e621ea65b7b2..6e9e21da7398 100644
--- a/Common/include/toolboxes/geometry_toolbox.hpp
+++ b/Common/include/toolboxes/geometry_toolbox.hpp
@@ -206,6 +206,20 @@ inline void RotationMatrix(Scalar psi, Matrix& mat) {
   mat[1][1] = cosPsi;
 }
 
+/*! \return Principal rotation angle in radians, using only the z angle in 2D. */
+template <class Scalar>
+inline Scalar RotationAngle(unsigned short nDim, const Scalar* angles) {
+  Scalar cosAngle;
+  if (nDim == 2) {
+    cosAngle = cos(angles[2]);
+  } else {
+    Scalar rotation[3][3];
+    RotationMatrix(angles[0], angles[1], angles[2], rotation);
+    cosAngle = 0.5 * (rotation[0][0] + rotation[1][1] + rotation[2][2] - 1.0);
+  }
+  return acos(std::max(Scalar(-1.0), std::min(Scalar(1.0), cosAngle)));
+}
+
 /*! \brief Apply a rotation matrix (R) about origin (O) to a point at
  *         distance (d) from it to obtain new coordinate (c). */
 template <class Scalar, int nDim>
diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp
index acafc2b7c025..446b98691677 100644
--- a/Common/src/geometry/CGeometry.cpp
+++ b/Common/src/geometry/CGeometry.cpp
@@ -127,9 +127,6 @@ CGeometry::~CGeometry() {
   delete[] bufD_PeriodicRecv;
   delete[] bufD_PeriodicSend;
 
-  delete[] bufS_PeriodicRecv;
-  delete[] bufS_PeriodicSend;
-
   delete[] req_PeriodicSend;
   delete[] req_PeriodicRecv;
 
@@ -1155,7 +1152,73 @@ void CGeometry::MatchPeriodic(const CConfig* config, unsigned short val_periodic
   delete[] Buffer_Recv_Marker;
 }
 
+void CGeometry::PreprocessPeriodicGeometry(const CConfig& config) {
+  const auto nPairs = config.GetnMarker_Periodic() / 2;
+  if (nPairs == 0) {
+    periodicPairIntersections = false;
+    periodicEdgeWeight.clear();
+    return;
+  }
+
+  vector<unsigned short> firstPair(nPoint, 0);
+  vector<bool> rotating(nPoint, false);
+  unsigned short local[2] = {}, global[2] = {};
+  for (auto iMarker = 0u; iMarker < nMarker; ++iMarker) {
+    if (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) continue;
+    const auto pair = (config.GetMarker_All_PerBound(iMarker) - 1) % nPairs + 1;
+    const auto rotation =
+        GeometryToolbox::HasRotation(config.GetPeriodicRotAngles(config.GetMarker_All_TagBound(iMarker)));
+    for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; ++iVertex) {
+      const auto iPoint = vertex[iMarker][iVertex]->GetNode();
+      if (nodes->GetVertex(iPoint, iMarker) < 0) continue;
+      if (firstPair[iPoint] != 0 && firstPair[iPoint] != pair) {
+        local[0] = 1;
+        if (rotation || rotating[iPoint]) local[1] = 1;
+      }
+      if (firstPair[iPoint] == 0) firstPair[iPoint] = pair;
+      rotating[iPoint] = rotating[iPoint] || rotation;
+    }
+  }
+  SU2_MPI::Allreduce(local, global, 2, MPI_UNSIGNED_SHORT, MPI_MAX, SU2_MPI::GetComm());
+  periodicPairIntersections = global[0] != 0;
+  if (global[1])
+    SU2_MPI::Error("Intersecting periodic pairs involving rotation are not supported yet.", CURRENT_FUNCTION);
+
+  periodicEdgeWeight.assign(nEdge, 1.0);
+  for (auto iEdge = 0ul; iEdge < nEdge; ++iEdge) {
+    const auto iPoint = edges->GetNode(iEdge, 0), jPoint = edges->GetNode(iEdge, 1);
+    if (!nodes->GetPeriodicBoundary(iPoint) || !nodes->GetPeriodicBoundary(jPoint)) continue;
+    passivedouble weight = 1.0;
+    for (auto iMarker = 0u; iMarker < nMarker; ++iMarker) {
+      if (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY || nodes->GetVertex(iPoint, iMarker) < 0 ||
+          nodes->GetVertex(jPoint, iMarker) < 0)
+        continue;
+      const auto tag = config.GetMarker_All_TagBound(iMarker);
+      const auto donor = config.GetMarker_Periodic_Donor(tag);
+      if (donor < nMarker && nodes->GetVertex(iPoint, donor) >= 0 && nodes->GetVertex(jPoint, donor) >= 0) {
+        /*--- An edge along the rotation axis is shared by all N sectors,
+         * rather than by two copies of each of the two periodic faces. ---*/
+        if (config.GetMarker_All_PerBound(iMarker) > config.GetnMarker_Periodic() / 2) continue;
+        const auto* angles = config.GetPeriodicRotAngles(tag);
+        const auto angle = SU2_TYPE::GetValue(GeometryToolbox::RotationAngle(nDim, angles));
+        constexpr passivedouble angleTol = 1e-4;
+        const auto nCopy = angle > angleTol ? static_cast<int>(2 * PI_NUMBER / angle + 0.5) : 0;
+        if (fabs(nCopy * angle - 2 * PI_NUMBER) > angleTol)
+          SU2_MPI::Error(
+              "Periodic points on the rotation axis require an angle that divides 360 degrees (marker " + tag + ").",
+              CURRENT_FUNCTION);
+        weight *= 1.0 / nCopy;
+      } else {
+        weight *= 0.5;
+      }
+    }
+    periodicEdgeWeight[iEdge] = weight;
+  }
+}
+
 void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) {
+  geometry->PreprocessPeriodicGeometry(*config);
+
   /*--- We start with the send and receive lists already available in
    the form of stored periodic point-donor pairs. We will loop through
    these markers and establish the neighboring ranks and number of
@@ -1313,9 +1376,6 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) {
   bufD_PeriodicSend = nullptr;
   bufD_PeriodicRecv = nullptr;
 
-  bufS_PeriodicSend = nullptr;
-  bufS_PeriodicRecv = nullptr;
-
   /*--- Allocate memory for the MPI requests if we need to communicate. ---*/
 
   if (nPeriodicSend > 0) {
@@ -1518,12 +1578,6 @@ void CGeometry::AllocatePeriodicComms(unsigned short countPerPeriodicPoint) {
 
     delete[] bufD_PeriodicRecv;
     bufD_PeriodicRecv = new su2double[nRecv]();
-
-    delete[] bufS_PeriodicSend;
-    bufS_PeriodicSend = new unsigned short[nSend]();
-
-    delete[] bufS_PeriodicRecv;
-    bufS_PeriodicRecv = new unsigned short[nRecv]();
   }
   END_SU2_OMP_SAFE_GLOBAL_ACCESS
 }
@@ -1565,10 +1619,6 @@ void CGeometry::PostPeriodicRecvs(CGeometry* geometry, const CConfig* config, CO
         SU2_MPI::Irecv(&(static_cast<su2double*>(bufD_PeriodicRecv)[offset]), count, MPI_DOUBLE, source, tag,
                        SU2_MPI::GetComm(), &(req_PeriodicRecv[iRecv]));
         break;
-      case COMM_TYPE::UNSIGNED_SHORT:
-        SU2_MPI::Irecv(&(static_cast<unsigned short*>(bufS_PeriodicRecv)[offset]), count, MPI_UNSIGNED_SHORT, source,
-                       tag, SU2_MPI::GetComm(), &(req_PeriodicRecv[iRecv]));
-        break;
       default:
         SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION);
         break;
@@ -1612,10 +1662,6 @@ void CGeometry::PostPeriodicSends(CGeometry* geometry, const CConfig* config, CO
         SU2_MPI::Isend(&(static_cast<su2double*>(bufD_PeriodicSend)[offset]), count, MPI_DOUBLE, dest, tag,
                        SU2_MPI::GetComm(), &(req_PeriodicSend[val_iSend]));
         break;
-      case COMM_TYPE::UNSIGNED_SHORT:
-        SU2_MPI::Isend(&(static_cast<unsigned short*>(bufS_PeriodicSend)[offset]), count, MPI_UNSIGNED_SHORT, dest, tag,
-                       SU2_MPI::GetComm(), &(req_PeriodicSend[val_iSend]));
-        break;
       default:
         SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION);
         break;
@@ -1636,9 +1682,6 @@ void CGeometry::PostPeriodicSends(CGeometry* geometry, const CConfig* config, CO
     case COMM_TYPE::DOUBLE:
       parallelCopy(myFinal - myStart, &bufD_PeriodicSend[myStart], &bufD_PeriodicRecv[iRecv]);
       break;
-    case COMM_TYPE::UNSIGNED_SHORT:
-      parallelCopy(myFinal - myStart, &bufS_PeriodicSend[myStart], &bufS_PeriodicRecv[iRecv]);
-      break;
     default:
       SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION);
       break;
@@ -1648,6 +1691,7 @@ void CGeometry::PostPeriodicSends(CGeometry* geometry, const CConfig* config, CO
 }
 
 void CGeometry::SetEdges() {
+  periodicEdgeWeight.clear();
   nEdge = 0;
   for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
     for (auto iNode = 0u; iNode < nodes->GetnPoint(iPoint); iNode++) {
diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp
index d3bca729bf5b..e02d613962ec 100644
--- a/Common/src/geometry/CMultiGridGeometry.cpp
+++ b/Common/src/geometry/CMultiGridGeometry.cpp
@@ -1157,7 +1157,14 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co
       auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren);
       if (fine_grid->nodes->GetBoundary(iFinePoint)) {
         nodes->SetBoundary(iCoarsePoint, nMarker);
-        break;
+        nodes->SetPeriodicBoundary(iCoarsePoint, nodes->GetPeriodicBoundary(iCoarsePoint) ||
+                                                     fine_grid->nodes->GetPeriodicBoundary(iFinePoint));
+        nodes->SetPhysicalBoundary(iCoarsePoint, nodes->GetPhysicalBoundary(iCoarsePoint) ||
+                                                     fine_grid->nodes->GetPhysicalBoundary(iFinePoint));
+        nodes->SetSolidBoundary(
+            iCoarsePoint, nodes->GetSolidBoundary(iCoarsePoint) || fine_grid->nodes->GetSolidBoundary(iFinePoint));
+        nodes->SetViscousBoundary(
+            iCoarsePoint, nodes->GetViscousBoundary(iCoarsePoint) || fine_grid->nodes->GetViscousBoundary(iFinePoint));
       }
     }
 
diff --git a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp
index 296f433ff40f..7df44db54b9a 100644
--- a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp
+++ b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp
@@ -32,6 +32,32 @@
 
 namespace detail {
 
+/*! \brief Visit the upper triangle and the extra xz entry used by the 3D LS factorization. */
+template<size_t nDim, class Function>
+FORCEINLINE void forEachLeastSquaresNormalEntry(Function function) {
+  for (size_t iDim = 0; iDim < nDim; ++iDim)
+    for (size_t jDim = iDim; jDim < nDim; ++jDim)
+      function(iDim, jDim, iDim, jDim);
+  if (nDim == 3) function(2, 1, 0, 2);
+}
+
+/*! \brief Rotate the symmetric normal matrix and pack it in the gradient storage convention. */
+template<size_t nDim, class MatrixType>
+FORCEINLINE void rotateLeastSquaresNormalMatrix(size_t iPoint, const MatrixType& matrix,
+                                               const su2double (&rotation)[nDim][nDim], su2double* packed) {
+  su2double rotated[nDim][nDim] = {};
+  for (size_t iDim = 0; iDim < nDim; ++iDim)
+    for (size_t jDim = iDim; jDim < nDim; ++jDim)
+      for (size_t kDim = 0; kDim < nDim; ++kDim)
+        for (size_t lDim = 0; lDim < nDim; ++lDim)
+          rotated[iDim][jDim] += rotation[iDim][kDim] * matrix(iPoint, min(kDim,lDim), max(kDim,lDim)) *
+                                rotation[jDim][lDim];
+  std::fill_n(packed, nDim * nDim, su2double(0));
+  forEachLeastSquaresNormalEntry<nDim>([&](size_t iDim, size_t jDim, size_t row, size_t col) {
+    packed[iDim * nDim + jDim] = rotated[row][col];
+  });
+}
+
 /*!
  * \brief Prepare Smatrix for 2D.
  * \ingroup FvmAlgos
@@ -231,8 +257,9 @@ void computeGradientsLeastSquares(CSolver* solver,
         Rmatrix(iPoint, iDim, jDim) = 0.0;
 
 
-    for (auto jPoint : nodes->GetPoints(iPoint))
-    {
+    for (auto iNeigh = 0u; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh) {
+      const auto jPoint = nodes->GetPoint(iPoint, iNeigh);
+      const auto iEdge = nodes->GetEdge(iPoint, iNeigh);
       const auto coord_j = geometry.nodes->GetCoord(jPoint);
       AD::SetPreaccIn(coord_j, nDim);
 
@@ -253,13 +280,12 @@ void computeGradientsLeastSquares(CSolver* solver,
       if (weight > 0.0)
       {
         weight = 1.0 / weight;
+        if (geometry.HasPeriodicPairIntersections())
+          weight *= geometry.GetPeriodicEdgeWeight(iEdge);
 
-        for (size_t iDim = 0; iDim < nDim; ++iDim)
-          for (size_t jDim = iDim; jDim < nDim; ++jDim)
-            Rmatrix(iPoint,iDim,jDim) += dist_ij[iDim]*dist_ij[jDim]*weight;
-
-        if (nDim == 3)
-          Rmatrix(iPoint,2,1) += dist_ij[0]*dist_ij[nDim-1]*weight;
+        forEachLeastSquaresNormalEntry<nDim>([&](size_t iDim, size_t jDim, size_t row, size_t col) {
+          Rmatrix(iPoint,iDim,jDim) += dist_ij[row]*dist_ij[col]*weight;
+        });
 
         /*--- Entries of c:= transpose(A)*b ---*/
 
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
index ae07bef84552..73ca86c3da83 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
@@ -791,8 +791,9 @@ class CFVMFlowSolverBase : public CSolver {
       jPoint_UndLapl[iPoint] = msw ? 1 : 0;
 
       /*--- Loop over the neighbors of point i. ---*/
-      for (auto jPoint : geometry->nodes->GetPoints(iPoint))
-      {
+      for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+        const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+        const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
         bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint);
 
         /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/
@@ -807,8 +808,10 @@ class CFVMFlowSolverBase : public CSolver {
           iPoint_UndLapl[iPoint] = fmax(iPoint_UndLapl[iPoint], fabs(sensVar_j - sensVar_i) / fmin(sensVar_j, sensVar_i));
         } else {
           /*--- Jameson dissipation sensor, add variable difference and variable sum. ---*/
-          iPoint_UndLapl[iPoint] += sensVar_j - sensVar_i;
-          jPoint_UndLapl[iPoint] += sensVar_j + sensVar_i;
+          const auto weight = geometry->HasPeriodicPairIntersections() ?
+                                  geometry->GetPeriodicEdgeWeight(iEdge) : 1.0;
+          iPoint_UndLapl[iPoint] += weight * (sensVar_j - sensVar_i);
+          jPoint_UndLapl[iPoint] += weight * (sensVar_j + sensVar_i);
         }
       }
 
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 3a46fc8c4db6..f0e4adf868ac 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -260,6 +260,13 @@ void CFVMFlowSolverBase<V, R>::CommunicateInitialState(CGeometry* geometry, cons
     InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
     CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
   }
+
+  /*--- The periodic communication updates the number of neighbors of the owned points only, update the halos. ---*/
+
+  if (config->GetnMarker_Periodic() > 0) {
+    geometry->InitiateComms(geometry, config, MPI_QUANTITIES::NEIGHBORS);
+    geometry->CompleteComms(geometry, config, MPI_QUANTITIES::NEIGHBORS);
+  }
   SetImplicitPeriodic(euler_implicit);
   SetRotatePeriodic(true);
 
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
index e35314b842ca..af50e7bb4997 100644
--- a/SU2_CFD/include/solvers/CSolver.hpp
+++ b/SU2_CFD/include/solvers/CSolver.hpp
@@ -134,6 +134,8 @@ class CSolver {
 
   /*--- End variables that need to go. ---*/
 
+  su2activevector periodicNeighborCount; /*!< \brief Partial periodic stencil count, including shared edges. */
+
   su2activevector iPoint_UndLapl;  /*!< \brief Auxiliary variable for the undivided Laplacians. */
   su2activevector jPoint_UndLapl;  /*!< \brief Auxiliary variable for the undivided Laplacians. */
 
diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
index b1c61ad377de..b27825ebac8b 100644
--- a/SU2_CFD/src/solvers/CEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -2397,20 +2397,25 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *co
       nodes->SetUnd_Lapl(iPoint, iVar, 0.0);
 
     /*--- Loop over the neighbors of point i. ---*/
-    for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+    for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+      const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+      const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
 
       bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint);
 
       /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/
       if (boundary_i && !boundary_j) continue;
 
+      const auto weight = geometry->HasPeriodicPairIntersections() ?
+                              geometry->GetPeriodicEdgeWeight(iEdge) : 1.0;
+
       /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/
 
       for (unsigned short iVar = 0; iVar < nVar; iVar++)
-        nodes->AddUnd_Lapl(iPoint, iVar, nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar));
+        nodes->AddUnd_Lapl(iPoint, iVar, weight * (nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar)));
 
       su2double Pressure_j = nodes->GetPressure(jPoint);
-      nodes->AddUnd_Lapl(iPoint, nVar-1, Pressure_j-Pressure_i);
+      nodes->AddUnd_Lapl(iPoint, nVar-1, weight * (Pressure_j-Pressure_i));
     }
   }
   END_SU2_OMP_FOR
diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp
index e23d25af07c2..2aacbf2cc4b3 100644
--- a/SU2_CFD/src/solvers/CHeatSolver.cpp
+++ b/SU2_CFD/src/solvers/CHeatSolver.cpp
@@ -148,12 +148,16 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, const CSolver* fl
     ghostNodes = make_unique<CHeatVariable>(Solution_Inf[0], maxMarkerVertices, nDim, nVar, config);
   }
 
-  /*--- Communicate and store volume and the number of neighbors for any dual CVs that lie on on periodic markers. ---*/
-  for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) {
-    InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME);
-    CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME);
-    InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
-    CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
+  /*--- Communicate and store volume and the number of neighbors for any dual CVs that lie on on periodic markers.
+   * With a flow solver on the same geometry this was already done by the flow solver, and the values are
+   * accumulated, so it must not be done twice. ---*/
+  if (!flow) {
+    for (auto iPeriodic = 1u; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) {
+      InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME);
+      CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME);
+      InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
+      CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
+    }
   }
   /*--- Store if implicit scheme is used. This has implications on the Residual and Jacobian handling for periodic
    * boundaries  ---*/
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index 407e1c9bc7a3..90353cbe9745 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -207,7 +207,7 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config,
       break;
     case PERIODIC_NEIGHBORS:
       COUNT_PER_POINT  = 1;
-      MPI_TYPE         = COMM_TYPE::UNSIGNED_SHORT;
+      MPI_TYPE         = COMM_TYPE::DOUBLE;
       break;
     case PERIODIC_RESIDUAL:
       COUNT_PER_POINT  = nVar + nVar*nVar + 1;
@@ -379,6 +379,20 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
     SU2_MPI::Error("The NEMO solvers do not support rotational periodicity yet.", CURRENT_FUNCTION);
   }
 
+  const bool intersectingPairs = geometry->HasPeriodicPairIntersections();
+  if (commType == PERIODIC_NEIGHBORS && intersectingPairs && val_periodic_index == 1) {
+    SU2_OMP_SAFE_GLOBAL_ACCESS(periodicNeighborCount.resize(geometry->GetnPoint());)
+    SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
+    for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) {
+      periodicNeighborCount(iPoint) = 0;
+      for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+        const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
+        periodicNeighborCount(iPoint) += geometry->GetPeriodicEdgeWeight(iEdge);
+      }
+    }
+    END_SU2_OMP_FOR
+  }
+
   /*--- Local variables ---*/
 
   bool boundary_i, boundary_j;
@@ -447,7 +461,6 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
   su2double *bufDSend = geometry->bufD_PeriodicSend;
 
-  unsigned short *bufSSend = geometry->bufS_PeriodicSend;
 
   /*--- Handle the different types of gradient and limiter. ---*/
 
@@ -570,20 +583,26 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
           case PERIODIC_NEIGHBORS:
 
+            if (intersectingPairs) {
+              bufDSend[buf_offset] = periodicNeighborCount(iPoint);
+              break;
+            }
+
             nNeighbor = 0;
-            for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+            for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+              const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
 
               /*--- Check if this neighbor lies on the periodic face so
                that we avoid double counting neighbors on both sides. If
                not, increment the count of neighbors for the donor. ---*/
 
-              if (!geometry->nodes->GetPeriodicBoundary(jPoint))
+              if (geometry->GetPeriodicEdgeWeight(iEdge) == 1.0)
                 nNeighbor++;
             }
 
             /*--- Store the number of neighbors in bufffer. ---*/
 
-            bufSSend[buf_offset] = nNeighbor;
+            bufDSend[buf_offset] = nNeighbor;
 
             break;
 
@@ -683,6 +702,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
           case PERIODIC_LAPLACIAN:
 
+            if (intersectingPairs) {
+              const auto* undLapl = base_nodes->GetUndivided_Laplacian(iPoint);
+              for (auto iField = 0u; iField < nVar; ++iField)
+                bufDSend[buf_offset + iField] = undLapl[iField];
+              if (rotate_periodic) Rotate(zeros, &undLapl[1], &bufDSend[buf_offset + 1]);
+              break;
+            }
+
             /*--- For JST, the undivided Laplacian must be computed
              consistently by using the complete control volume info
              from both sides of the periodic face. ---*/
@@ -690,12 +717,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
             for (iVar = 0; iVar < nVar; iVar++)
               Und_Lapl[iVar] = 0.0;
 
-            for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+            for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+              const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+              const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
 
               /*--- Avoid periodic boundary points so that we do not
                duplicate edges on both sides of the periodic BC. ---*/
 
-              if (!geometry->nodes->GetPeriodicBoundary(jPoint)) {
+              if (geometry->GetPeriodicEdgeWeight(iEdge) == 1.0) {
 
                 /*--- Solution differences ---*/
 
@@ -749,6 +778,11 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
             break;
 
           case PERIODIC_SENSOR: {
+            if (intersectingPairs) {
+              bufDSend[buf_offset] = iPoint_UndLapl(iPoint);
+              bufDSend[buf_offset + 1] = jPoint_UndLapl(iPoint);
+              break;
+            }
             const bool msw = config->GetKind_Upwind_Flow() == UPWIND::MSW;
 
             /*--- For the centered schemes, the sensor must be computed
@@ -756,12 +790,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
              on both sides of the periodic face. ---*/
 
             Sensor_i = 0; Sensor_j = 0;
-            for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+            for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+              const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+              const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
 
               /*--- Avoid halos and boundary points so that we don't
                duplicate edges on both sides of the periodic BC. ---*/
 
-              if (geometry->nodes->GetPeriodicBoundary(jPoint)) continue;
+              if (geometry->GetPeriodicEdgeWeight(iEdge) < 1.0) continue;
 
               /*--- Use density instead of pressure for incomp. flows. ---*/
 
@@ -781,7 +817,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
               if ((!boundary_i || boundary_j) && geometry->nodes->GetDomain(iPoint)) {
                 if (msw) {
-                  Sensor_i = fmax(Sensor_i, fabs(Pressure_j - Pressure_i)) / fmin(Pressure_i, Pressure_j);
+                  Sensor_i = fmax(Sensor_i, fabs(Pressure_j - Pressure_i) / fmin(Pressure_i, Pressure_j));
                 } else {
                   Sensor_i += (Pressure_j - Pressure_i);
                   Sensor_j += (Pressure_i + Pressure_j);
@@ -859,6 +895,28 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
             /*--- Set a flag for unweighted or weighted least-squares. ---*/
 
+            if (intersectingPairs) {
+              if (nDim == 2)
+                detail::rotateLeastSquaresNormalMatrix(iPoint, base_nodes->GetRmatrix(), rotMatrix2D, &bufDSend[buf_offset]);
+              else
+                detail::rotateLeastSquaresNormalMatrix(iPoint, base_nodes->GetRmatrix(), rotMatrix3D, &bufDSend[buf_offset]);
+              buf_offset += static_cast<unsigned long>(nDim) * nDim;
+              for (auto iField = 0u; iField < ICOUNT; ++iField)
+                Rotate(zeros, gradient[iPoint][iField], rotBlock[iField]);
+              if (rotate_periodic) {
+                for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) {
+                  su2double velocity[3] = {}, rotatedVelocity[3] = {};
+                  for (auto jDim = 0u; jDim < nDim; ++jDim) velocity[jDim] = rotBlock[1+jDim][iCoordinate];
+                  Rotate(zeros, velocity, rotatedVelocity);
+                  for (auto jDim = 0u; jDim < nDim; ++jDim) rotBlock[1+jDim][iCoordinate] = rotatedVelocity[jDim];
+                }
+              }
+              for (auto iField = 0u; iField < ICOUNT; ++iField)
+                for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate)
+                  bufDSend[buf_offset++] = rotBlock[iField][iCoordinate];
+              break;
+            }
+
             switch(commType) {
               case PERIODIC_SOL_ULS:
               case PERIODIC_SOL_ULS_R:
@@ -899,12 +957,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
             r11 = 0.0;   r12 = 0.0;   r22 = 0.0;
             r13 = 0.0; r23_a = 0.0; r23_b = 0.0;  r33 = 0.0;
 
-            for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+            for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+              const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+              const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
 
               /*--- Avoid periodic boundary points so that we do not
                duplicate edges on both sides of the periodic BC. ---*/
 
-              if (!geometry->nodes->GetPeriodicBoundary(jPoint)) {
+              if (geometry->GetPeriodicEdgeWeight(iEdge) == 1.0) {
 
                 /*--- Get coordinates for the neighbor point. ---*/
 
@@ -1153,7 +1213,6 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry,
 
   const su2double *bufDRecv = geometry->bufD_PeriodicRecv;
 
-  const unsigned short *bufSRecv = geometry->bufS_PeriodicRecv;
 
   /*--- Handle the different types of gradient and limiter. ---*/
 
@@ -1230,13 +1289,13 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry,
               break;
 
             case PERIODIC_NEIGHBORS:
-
-              /*--- Store the extra neighbors on the periodic face. ---*/
-
-              nNeighbor = (geometry->nodes->GetnNeighbor(iPoint) +
-                           bufSRecv[buf_offset]);
-              geometry->nodes->SetnNeighbor(iPoint, nNeighbor);
-
+              if (geometry->HasPeriodicPairIntersections()) {
+                periodicNeighborCount(iPoint) += bufDRecv[buf_offset];
+                geometry->nodes->SetnNeighbor(iPoint, SU2_TYPE::Int(periodicNeighborCount(iPoint) + 0.5));
+              } else {
+                nNeighbor = geometry->nodes->GetnNeighbor(iPoint) + SU2_TYPE::Int(bufDRecv[buf_offset]);
+                geometry->nodes->SetnNeighbor(iPoint, nNeighbor);
+              }
               break;
 
             case PERIODIC_RESIDUAL:
@@ -2394,7 +2453,9 @@ void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config)
       base_nodes->SetUnd_Lapl(iPoint, iVar, 0.0);
 
     /*--- Loop over the neighbors of point i. ---*/
-    for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+    for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+      const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+      const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
 
       bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint);
 
@@ -2405,6 +2466,7 @@ void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config)
 
       for (unsigned short iVar = 0; iVar < nVar; iVar++) {
         su2double delta = base_nodes->GetSolution(jPoint,iVar)-base_nodes->GetSolution(iPoint,iVar);
+        if (geometry->HasPeriodicPairIntersections()) delta *= geometry->GetPeriodicEdgeWeight(iEdge);
         base_nodes->AddUnd_Lapl(iPoint, iVar, delta);
       }
     }
diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index 21c63694f567..d298966939fb 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -149,7 +149,7 @@ def main():
     poiseuille_profile.cfg_dir   = "navierstokes/poiseuille"
     poiseuille_profile.cfg_file  = "profile_poiseuille.cfg"
     poiseuille_profile.test_iter = 10
-    poiseuille_profile.test_vals         = [-12.004080, -7.638116, -0.000000, 2.089953]
+    poiseuille_profile.test_vals = [-12.004094, -7.640551, -0.000000, 2.089953]
     poiseuille_profile.test_vals_aarch64 = [-12.004276, -7.636719, -0.000000, 2.089953]
     test_list.append(poiseuille_profile)
 
@@ -178,7 +178,7 @@ def main():
     rae2822_sst.cfg_dir   = "rans/rae2822"
     rae2822_sst.cfg_file  = "turb_SST_RAE2822.cfg"
     rae2822_sst.test_iter = 20
-    rae2822_sst.test_vals = [-1.461774, 5.883885, 0.630267, 0.020404, 100.000000]
+    rae2822_sst.test_vals = [-1.519695, 5.888685, 0.614691, 0.020767, 100.000000]
     test_list.append(rae2822_sst)
 
     # RAE2822 SST_SUST
@@ -602,7 +602,7 @@ def main():
     Jones_tc_restart.cfg_dir   = "turbomachinery/APU_turbocharger"
     Jones_tc_restart.cfg_file  = "Jones_restart.cfg"
     Jones_tc_restart.test_iter = 5
-    Jones_tc_restart.test_vals = [-11.959252, -12.213981, -19.314222, -13.600248, -19.076384, -13.444335, 73286.000000, 73286.000000, 0.020056, 82.286000]
+    Jones_tc_restart.test_vals = [-11.961836, -12.211509, -19.315122, -13.598954, -19.081506, -13.442877, 73286.000000, 73286.000000, 0.020056, 82.286000]
     test_list.append(Jones_tc_restart)
 
     # 2D axial stage
@@ -711,7 +711,7 @@ def main():
     bars_SST_2D.cfg_dir   = "sliding_interface/bars_SST_2D"
     bars_SST_2D.cfg_file  = "bars.cfg"
     bars_SST_2D.test_iter = 13
-    bars_SST_2D.test_vals = [13.000000, -0.455981, -1.541069]
+    bars_SST_2D.test_vals = [13.000000, -0.455976, -1.539981]
     bars_SST_2D.multizone = True
     test_list.append(bars_SST_2D)
 
diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg
new file mode 100644
index 000000000000..6ab487b08981
--- /dev/null
+++ b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg
@@ -0,0 +1,79 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: 2D pin array, periodic flow driven by a body force, with   %
+%                   the weakly coupled heat equation (flow and heat solver on  %
+%                   the same mesh, the heat solver does not change the flow)   %
+% File Version 8.5.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+SOLVER= INC_NAVIER_STOKES
+KIND_TURB_MODEL= NONE
+RESTART_SOL= NO
+%
+INC_ENERGY_EQUATION= NO
+WEAKLY_COUPLED_HEAT_EQUATION= YES
+
+% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------%
+%
+INC_DENSITY_MODEL= CONSTANT
+INC_DENSITY_INIT= 1045.0
+INC_VELOCITY_INIT= ( 0.0001, 0.0, 0.0 )
+INC_TEMPERATURE_INIT= 338.0
+INC_NONDIM= DIMENSIONAL
+%
+SPECIFIC_HEAT_CP= 3540.0
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+MU_CONSTANT= 0.001385
+CONDUCTIVITY_MODEL= CONSTANT_PRANDTL
+PRANDTL_LAM= 11.7
+%
+BODY_FORCE= YES
+BODY_FORCE_VECTOR= ( 0.5, 0.0, 0.0 )
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_SYM= ( fluid_symmetry )
+MARKER_PERIODIC= ( fluid_inlet, fluid_outlet, 0.0,0.0,0.0, 0.0,0.0,0.0, 0.0111544,0.0,0.0 )
+MARKER_HEATFLUX= ( fluid_pin1_interface, 0.0, fluid_pin2_interface, 0.0, fluid_pin3_interface, 0.0 )
+%
+MARKER_MONITORING= ( fluid_pin2_interface )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= GREEN_GAUSS
+CFL_NUMBER= 1000
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1e-8
+LINEAR_SOLVER_ITER= 20
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= FDS
+MUSCL_FLOW= YES
+SLOPE_LIMITER_FLOW= NONE
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_HEAT= SCALAR_UPWIND
+TIME_DISCRE_HEAT= EULER_IMPLICIT
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+CONV_FIELD= RMS_PRESSURE
+CONV_RESIDUAL_MINVAL= -16
+CONV_STARTITER= 10
+ITER= 2000
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FORMAT= SU2
+MESH_FILENAME= fluid.su2
+OUTPUT_WRT_FREQ= 9999
+SCREEN_OUTPUT= ( INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, DRAG )
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 30824edc691b..a0a2eeda708b 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -425,7 +425,7 @@ def main():
     poiseuille_profile.cfg_dir   = "navierstokes/poiseuille"
     poiseuille_profile.cfg_file  = "profile_poiseuille.cfg"
     poiseuille_profile.test_iter = 10
-    poiseuille_profile.test_vals = [-12.003186, -7.623957, -0.000000, 2.089953]
+    poiseuille_profile.test_vals = [-12.003181, -7.626400, -0.000000, 2.089953]
     poiseuille_profile.test_vals_aarch64 = [-12.007498, -7.226926, -0.000000, 2.089953]
     poiseuille_profile.tol       = [0.001, 0.001, 1e-5, 1e-5, 1e-5]
     test_list.append(poiseuille_profile)
@@ -846,6 +846,14 @@ def main():
     inc_heatTransfer_BC.test_vals = [-8.904266, -7.745636, -8.064003, -1.079197, -1671.000000]
     test_list.append(inc_heatTransfer_BC)
 
+    # 2D pin array, periodic with a body force, flow and weakly coupled heat equation
+    inc_periodic_weak_heat = TestCase('inc_periodic_weak_heat')
+    inc_periodic_weak_heat.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d"
+    inc_periodic_weak_heat.cfg_file = "periodic_weak_heat.cfg"
+    inc_periodic_weak_heat.test_iter = 10
+    inc_periodic_weak_heat.test_vals = [-4.712559, -5.747243, -5.751727, 873.702714]
+    test_list.append(inc_periodic_weak_heat)
+
     ############################
     ### Incompressible RANS  ###
     ############################
@@ -957,7 +965,7 @@ def main():
     turbmod_sa_neg_rae2822.cfg_dir   = "turbulence_models/sa/rae2822"
     turbmod_sa_neg_rae2822.cfg_file  = "turb_SA_NEG_RAE2822.cfg"
     turbmod_sa_neg_rae2822.test_iter = 10
-    turbmod_sa_neg_rae2822.test_vals = [1.275994, 1.062588, -2.067822, 1.300662, 0.480881, 0.000000]
+    turbmod_sa_neg_rae2822.test_vals = [1.275906, 1.062762, -2.067335, 1.301619, 0.481422, 0.000000]
     turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.345593, 1.448310, 1.208721, -0.846597, 1.248410, 0.489117, 0.000000]
     test_list.append(turbmod_sa_neg_rae2822)
 
@@ -1243,7 +1251,7 @@ def main():
     sbs_backward_step.cfg_dir    = "backscatter/backward_step"
     sbs_backward_step.cfg_file   = "backwardStep.cfg"
     sbs_backward_step.test_iter  = 3
-    sbs_backward_step.test_vals  = [-6.352884, -3.465372, -5.507901, -3.906545, -9.506305, -6.365236, -6.331021, -6.331028]
+    sbs_backward_step.test_vals = [-6.352885, -3.465372, -5.507901, -3.906544, -9.506305, -6.365236, -6.331021, -6.331028]
     sbs_backward_step.unsteady   = True
     sbs_backward_step.decompress = True
     sbs_backward_step.grid_file  = "backward_step.su2"
@@ -1323,7 +1331,7 @@ def main():
     Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg"
     Aachen_3D_restart.test_iter = 5
     Aachen_3D_restart.tol = 0.00001
-    Aachen_3D_restart.test_vals = [-7.688483, -8.466346, -6.035067, -6.425967, -5.820125, -4.597817, -5.528297, -5.301466, -3.819862, -5.242356, -5.738805, -3.630202, -2.213240, -2.857919, -0.576864]
+    Aachen_3D_restart.test_vals = [-7.688401, -8.466746, -6.034652, -6.427491, -5.818347, -4.597687, -5.530255, -5.303027, -3.819060, -5.242491, -5.755143, -3.629921, -2.213103, -2.859285, -0.576337]
     test_list.append(Aachen_3D_restart)
 
     # Jones APU Turbocharger restart
@@ -1331,7 +1339,7 @@ def main():
     Jones_tc_restart.cfg_dir   = "turbomachinery/APU_turbocharger"
     Jones_tc_restart.cfg_file  = "Jones_restart.cfg"
     Jones_tc_restart.test_iter = 5
-    Jones_tc_restart.test_vals = [-11.905382, -12.212100, -19.212143, -13.508268, -19.078783, -13.446389, 73286.000000, 73286.000000, 0.020056, 82.286000]
+    Jones_tc_restart.test_vals = [-11.903628, -12.211875, -19.212536, -13.508088, -19.079442, -13.440735, 73286.000000, 73286.000000, 0.020056, 82.286000]
     test_list.append(Jones_tc_restart)
 
     # 2D axial stage
@@ -1439,7 +1447,7 @@ def main():
     bars_SST_2D.cfg_dir   = "sliding_interface/bars_SST_2D"
     bars_SST_2D.cfg_file  = "bars.cfg"
     bars_SST_2D.test_iter = 13
-    bars_SST_2D.test_vals = [13.000000, -0.457912, -1.541047]
+    bars_SST_2D.test_vals = [13.000000, -0.457907, -1.539958]
     bars_SST_2D.multizone = True
     test_list.append(bars_SST_2D)
 
@@ -1765,7 +1773,7 @@ def main():
     mms_fvm_vortex.cfg_dir   = "mms/fvm_euler"
     mms_fvm_vortex.cfg_file  = "inv_mms_vortex.cfg"
     mms_fvm_vortex.test_iter = 10
-    mms_fvm_vortex.test_vals = [-5.704300, -4.848072, 0.000000, 0.000000]
+    mms_fvm_vortex.test_vals = [-5.704011, -4.847997, 0.000000, 0.000000]
     mms_fvm_vortex.unsteady  = True
     test_list.append(mms_fvm_vortex)
 
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 250f96e66c12..ef8999cadd0e 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -315,7 +315,7 @@ def main():
     discadj_fsi2.cfg_dir   = "disc_adj_fsi/Airfoil_2d"
     discadj_fsi2.cfg_file  = "config.cfg"
     discadj_fsi2.test_iter = 8
-    discadj_fsi2.test_vals         = [-3.824633, 1.979516, -3.863368, 0.295450, 3.839800]
+    discadj_fsi2.test_vals = [-3.824660, 1.979488, -3.863368, 0.295450, 3.839800]
     discadj_fsi2.test_vals_aarch64 = [-3.824870, 1.979160, -3.863368, 0.295450, 3.839800]
     discadj_fsi2.tol       = 0.00001
     test_list.append(discadj_fsi2)
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index 1edd0378ab2a..e65b6586eae4 100755
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -216,7 +216,7 @@ def main():
     poiseuille_profile.cfg_dir   = "navierstokes/poiseuille"
     poiseuille_profile.cfg_file  = "profile_poiseuille.cfg"
     poiseuille_profile.test_iter = 10
-    poiseuille_profile.test_vals         = [-12.003115, -7.626023, -0.000000, 2.089953]
+    poiseuille_profile.test_vals = [-12.003109, -7.627917, -0.000000, 2.089953]
     poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262299, -0.000000, 2.089953] #last 4 columns
     test_list.append(poiseuille_profile)
 
@@ -957,7 +957,7 @@ def main():
     Aachen_3D_restart.cfg_dir   = "turbomachinery/Aachen_turbine"
     Aachen_3D_restart.cfg_file  = "aachen_3D_MP_restart.cfg"
     Aachen_3D_restart.test_iter = 5
-    Aachen_3D_restart.test_vals = [-7.689117, -8.476204, -6.036107, -6.427402, -5.831576, -4.599024, -5.530512, -5.305687, -3.821789, -5.243255, -5.754831, -3.631603, -2.213934, -2.895728, -0.578545]
+    Aachen_3D_restart.test_vals = [-7.689029, -8.475341, -6.035704, -6.426165, -5.830215, -4.598893, -5.529338, -5.301964, -3.820993, -5.243202, -5.732910, -3.631319, -2.213792, -2.876325, -0.578057]
     Aachen_3D_restart.enabled_with_asan = False
     test_list.append(Aachen_3D_restart)
 
@@ -966,7 +966,7 @@ def main():
     Jones_tc_restart.cfg_dir   = "turbomachinery/APU_turbocharger"
     Jones_tc_restart.cfg_file  = "Jones_restart.cfg"
     Jones_tc_restart.test_iter = 5
-    Jones_tc_restart.test_vals = [-11.905642, -12.212181, -19.212554, -13.508053, -19.069991, -13.439494, 73286.000000, 73286.000000, 0.020056, 82.286000]
+    Jones_tc_restart.test_vals = [-11.935390, -12.212295, -19.248052, -13.539427, -19.086000, -13.455882, 73286.000000, 73286.000000, 0.020056, 82.286000]
     test_list.append(Jones_tc_restart)
 
     # 2D axial stage
@@ -1086,7 +1086,7 @@ def main():
     bars_SST_2D.cfg_dir   = "sliding_interface/bars_SST_2D"
     bars_SST_2D.cfg_file  = "bars.cfg"
     bars_SST_2D.test_iter = 13
-    bars_SST_2D.test_vals = [13.000000, -0.456143, -1.541051]
+    bars_SST_2D.test_vals = [13.000000, -0.456139, -1.539962]
     bars_SST_2D.multizone = True
     test_list.append(bars_SST_2D)
 
@@ -1238,7 +1238,7 @@ def main():
     mms_fvm_vortex.cfg_dir   = "mms/fvm_euler"
     mms_fvm_vortex.cfg_file  = "inv_mms_vortex.cfg"
     mms_fvm_vortex.test_iter = 10
-    mms_fvm_vortex.test_vals = [-5.704300, -4.848072, 0.000000, 0.000000]
+    mms_fvm_vortex.test_vals = [-5.704011, -4.847997, 0.000000, 0.000000]
     mms_fvm_vortex.unsteady  = True
     test_list.append(mms_fvm_vortex)
 
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 288f8a47efd3..44bfc4d0073d 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -200,7 +200,7 @@ def main():
     discadj_trans_stator.cfg_dir   = "disc_adj_turbomachinery/transonic_stator_2D"
     discadj_trans_stator.cfg_file  = "transonic_stator.cfg"
     discadj_trans_stator.test_iter = 79
-    discadj_trans_stator.test_vals = [79.000000, -7.308168, -9.891121, -10.038669, -13.368502]
+    discadj_trans_stator.test_vals = [79.000000, -7.308179, -9.891140, -10.038680, -13.368512]
     discadj_trans_stator.test_vals_aarch64 = [79.000000, -7.308167, -9.891117, -10.038669, -13.368501]
     test_list.append(discadj_trans_stator)
 
diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py
index 22e9e55fdb03..f08b73f21d26 100644
--- a/TestCases/tutorials.py
+++ b/TestCases/tutorials.py
@@ -314,7 +314,7 @@ def main():
     tutorial_unst_naca0012.cfg_dir       = "../Tutorials/compressible_flow/Unsteady_NACA0012"
     tutorial_unst_naca0012.cfg_file      = "unsteady_naca0012.cfg"
     tutorial_unst_naca0012.test_iter     = 520
-    tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.292186, 0.000000, 0.304255, 0.792893, 0.001907, 0.010888]
+    tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.296539, 0.000000, 0.308341, 0.797028, 0.003155, 0.015433]
     tutorial_unst_naca0012.test_vals_aarch64 = [520, 0, -5.292359, 0, 0.284720, 0.766329, 0.000954, 0.007565]
     tutorial_unst_naca0012.unsteady      = True
     test_list.append(tutorial_unst_naca0012)
diff --git a/UnitTests/SU2_CFD/check_periodic_intersection_rejection.py b/UnitTests/SU2_CFD/check_periodic_intersection_rejection.py
new file mode 100644
index 000000000000..e1ac8eb0f751
--- /dev/null
+++ b/UnitTests/SU2_CFD/check_periodic_intersection_rejection.py
@@ -0,0 +1,18 @@
+"""Run the MPI-aborting geometry check in a separate process."""
+
+import subprocess
+import sys
+
+result = subprocess.run(
+    sys.argv[1:] + ["Invalid rotational intersection"],
+    stdout=subprocess.PIPE,
+    stderr=subprocess.STDOUT,
+    text=True,
+    timeout=60,
+)
+print(result.stdout)
+assert result.returncode != 0, "Rotational intersection was accepted"
+assert (
+    "Intersecting periodic pairs involving rotation are not supported yet."
+    in result.stdout
+)
diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp
index 72baf525ca68..95b37db5cf9e 100644
--- a/UnitTests/SU2_CFD/gradients.cpp
+++ b/UnitTests/SU2_CFD/gradients.cpp
@@ -26,6 +26,8 @@
  */
 
 #include "catch.hpp"
+#include "../UnitQuadTestCase.hpp"
+#include "../../Common/include/geometry/CMultiGridGeometry.hpp"
 #include "../../Common/include/geometry/CPhysicalGeometry.hpp"
 #include "../../Common/include/containers/container_decorators.hpp"
 #include "../../SU2_CFD/include/solvers/CSolver.hpp"
@@ -154,3 +156,359 @@ TEST_CASE("GG", "[Gradients]") { testGreenGauss<LinearFunction>(); }
 TEST_CASE("LS", "[Gradients]") { testLeastSquares<LinearFunction>(false); }
 
 TEST_CASE("WLS", "[Gradients]") { testLeastSquares<LinearFunction>(true); }
+
+TEST_CASE("Intersecting periodic stencils", "[Periodic][Gradients]") {
+  const auto nPairs = GENERATE(2u, 3u);
+  const bool diagonal = GENERATE(false, true);
+  const auto scheme = GENERATE(0u, 1u, 2u);
+  const bool msw = scheme == 1, incompressible = scheme == 2;
+  UnitQuadTestCase field;
+  const auto start = field.config_options.find("MARKER_HEATFLUX=");
+  const auto end = field.config_options.find("VISCOSITY_MODEL=");
+  field.config_options.replace(start, end - start, nPairs == 2 ? "MARKER_CUSTOM= (z_plus,z_minus)\n" : "");
+  std::string periodic = "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0, y_minus,y_plus, 0,0,0, 0,0,0, 0,1,0";
+  if (nPairs == 3) periodic += ", z_minus,z_plus, 0,0,0, 0,0,0, 0,0,1";
+  field.AddOption(periodic + ")");
+  field.AddOption(msw ? "CONV_NUM_METHOD_FLOW= MSW" : "CONV_NUM_METHOD_FLOW= JST");
+  field.AddOption("MUSCL_FLOW= NO\nNUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES");
+  if (incompressible) {
+    field.SetOption("SOLVER= INC_NAVIER_STOKES");
+    field.AddOption("FLUID_MODEL= CONSTANT_DENSITY");
+    field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION");
+  }
+  field.InitConfig();
+  field.InitGeometry(true);
+  if (diagonal) {
+    /*--- A triangulated plane adds NE/SW edges. The SW donor at a periodic
+     * corner is reached through two successive pairs, rather than one. ---*/
+    std::vector<std::vector<unsigned long>> neighbors(field.geometry->GetnPoint());
+    for (auto i = 0ul; i < neighbors.size(); ++i) {
+      for (auto j : field.geometry->nodes->GetPoints(i)) neighbors[i].push_back(j);
+      const auto* x = field.geometry->nodes->GetCoord(i);
+      for (auto j = 0ul; j < neighbors.size(); ++j) {
+        const auto* y = field.geometry->nodes->GetCoord(j);
+        if (fabs(y[0] - x[0] - 0.25) < 1e-12 && fabs(y[1] - x[1] - 0.25) < 1e-12 && fabs(y[2] - x[2]) < 1e-12) {
+          neighbors[i].push_back(j);
+          neighbors[j].push_back(i);
+        }
+      }
+    }
+    field.geometry->nodes->SetPoints(neighbors);
+    delete field.geometry->edges;
+    field.geometry->SetEdges();
+    for (auto i = 0ul; i < neighbors.size(); ++i) field.geometry->nodes->SetnNeighbor(i, neighbors[i].size());
+  }
+  for (auto pair = 1u; pair <= nPairs; ++pair) field.geometry->MatchPeriodic(field.config.get(), pair);
+  field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get());
+  CHECK(field.geometry->HasPeriodicPairIntersections() == (field.config->GetnMarker_Periodic() > 2));
+  field.InitSolver();
+  auto* solver = field.solver[FLOW_SOL];
+  auto* nodes = solver->GetNodes();
+  auto value = [](su2double x, su2double y, su2double z) {
+    return 10 + x * (1 - x) + 2 * y * (1 - y) + 3 * z * (1 - z);
+  };
+  unsigned long target = field.geometry->GetnPoint();
+  for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    nodes->SetPressure(i, value(x[0], x[1], x[2]));
+    for (auto v = 0u; v < solver->GetnVar(); ++v)
+      nodes->SetSolution(
+          i, v,
+          v == 0 ? value(x[0], x[1], x[2])
+                 : (v + 1 == solver->GetnVar() ? (incompressible ? 300.0 : value(x[0], x[1], x[2]) / 0.4) : 0));
+    if (field.geometry->nodes->GetDomain(i) && x[0] == 0 && x[1] == 0 && x[2] == (nPairs == 3 ? 0 : 0.5)) target = i;
+  }
+  const auto z = nPairs == 3 ? 0.0 : 0.5;
+  const auto center = value(0, 0, z);
+  std::vector<su2double> stencil = {value(0.25, 0, z),     value(0.75, 0, z),
+                                    value(0, 0.25, z),     value(0, 0.75, z),
+                                    value(0, 0, z + 0.25), value(0, 0, nPairs == 3 ? 0.75 : z - 0.25)};
+  if (diagonal) {
+    stencil.push_back(value(0.25, 0.25, z));
+    stencil.push_back(value(0.75, 0.75, z));
+  }
+  su2double numerator = 0, denominator = 0, maximum = 0;
+  for (const auto neighbor : stencil) {
+    numerator += neighbor - center;
+    denominator += neighbor + center;
+    maximum = std::max(maximum, fabs(neighbor - center) / std::min(center, neighbor));
+  }
+  field.config->SetGlobalParam(field.config->GetKind_Solver(), RUNTIME_FLOW_SYS);
+  {
+    SU2_OMP_PARALLEL {
+      solver->Preprocessing(field.geometry.get(), field.solver, field.config.get(), 0, 0, RUNTIME_FLOW_SYS, false);
+    }
+    END_SU2_OMP_PARALLEL
+  }
+  unsigned long localTarget = target < field.geometry->GetnPoint(), globalTarget = 0;
+  SU2_MPI::Allreduce(&localTarget, &globalTarget, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  REQUIRE(globalTarget == 1);
+  INFO("pairs " << nPairs << ", diagonal " << diagonal << ", MSW " << msw << ", incompressible " << incompressible);
+  if (target < field.geometry->GetnPoint()) {
+    CHECK(field.geometry->nodes->GetnNeighbor(target) == stencil.size());
+    CHECK(nodes->GetSensor(target) ==
+          Approx(msw ? maximum : (incompressible ? 0.0 : fabs(numerator) / denominator)).margin(1e-12));
+    if (!msw) CHECK(nodes->GetUndivided_Laplacian(target, 0) == Approx(numerator).margin(1e-12));
+  }
+  /*--- Periodic scalar field: compare the assembled primitive LS normal
+   * equations with the complete stencil, including the composed diagonal. ---*/
+  solver->SetRotatePeriodic(false);
+  auto& gradient = nodes->GetGradient_Primitive();
+  auto& matrix = nodes->GetRmatrix();
+  for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) nodes->SetPrimitive(i, v, value(x[0], x[1], x[2]));
+  }
+  SU2_OMP_PARALLEL {
+    computeGradientsLeastSquares(solver, MPI_QUANTITIES::PRIMITIVE_GRADIENT, PERIODIC_PRIM_LS, *field.geometry,
+                                 *field.config, true, nodes->GetPrimitive(), 0, solver->GetnPrimVarGrad(), -1, gradient,
+                                 matrix);
+  }
+  END_SU2_OMP_PARALLEL
+  if (target < field.geometry->GetnPoint()) {
+    CHECK(gradient(target, 0, 0) == Approx(0).margin(1e-12));
+    CHECK(gradient(target, 0, 1) == Approx(0).margin(1e-12));
+    CHECK(gradient(target, 0, 2) == Approx(nPairs == 3 ? 0 : 3 * (1 - 2 * z)).margin(1e-12));
+  }
+}
+
+template <size_t nDim>
+void checkRotatedLeastSquaresNormalMatrix() {
+  C3DDoubleMatrix matrix(1, nDim, nDim, 0.0);
+  const su2double vectors[][3] = {{1, 2, 3}, {2, -1, 1}, {-1, 1, 2}};
+  su2double rotation[nDim][nDim] = {};
+  if constexpr (nDim == 2)
+    GeometryToolbox::RotationMatrix(su2double(PI_NUMBER / 2), rotation);
+  else
+    GeometryToolbox::RotationMatrix(su2double(0), su2double(0), su2double(PI_NUMBER / 2), rotation);
+  su2double expected[nDim][nDim] = {};
+  for (const auto& vector : vectors) {
+    detail::forEachLeastSquaresNormalEntry<nDim>(
+        [&](size_t i, size_t j, size_t row, size_t col) { matrix(0, i, j) += vector[row] * vector[col]; });
+    su2double rotated[nDim] = {}, zeros[nDim] = {};
+    GeometryToolbox::Rotate(rotation, zeros, vector, rotated);
+    for (auto i = 0u; i < nDim; ++i)
+      for (auto j = 0u; j < nDim; ++j) expected[i][j] += rotated[i] * rotated[j];
+  }
+  su2double packed[nDim * nDim];
+  detail::rotateLeastSquaresNormalMatrix(0, matrix, rotation, packed);
+  for (auto i = 0u; i < nDim; ++i)
+    for (auto j = 0u; j < nDim; ++j) {
+      const auto ref = i <= j ? expected[i][j] : (nDim == 3 && i == 2 && j == 1 ? expected[0][nDim - 1] : 0.0);
+      CHECK(packed[i * nDim + j] == Approx(ref).margin(1e-12));
+    }
+}
+
+TEST_CASE("Rotated least-squares normal matrix", "[Periodic][Gradients]") {
+  checkRotatedLeastSquaresNormalMatrix<2>();
+  checkRotatedLeastSquaresNormalMatrix<3>();
+}
+
+TEST_CASE("2D periodic linear gradient", "[Periodic][Gradients]") {
+  const bool weighted = GENERATE(false, true);
+  UnitQuadTestCase field;
+  field.SetOption("MESH_FORMAT= RECTANGLE");
+  field.SetOption("MARKER_CUSTOM= (z_plus,z_minus)");
+  field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)");
+  field.InitConfig();
+  field.InitGeometry(true);
+  field.geometry->MatchPeriodic(field.config.get(), 1);
+  field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get());
+  field.InitSolver();
+  auto* solver = field.solver[FLOW_SOL];
+  auto* nodes = solver->GetNodes();
+  for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i)
+    for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v)
+      nodes->SetPrimitive(i, v, 3 + field.geometry->nodes->GetCoord(i, 1));
+  solver->SetRotatePeriodic(false);
+  SU2_OMP_PARALLEL {
+    computeGradientsLeastSquares(solver, MPI_QUANTITIES::PRIMITIVE_GRADIENT,
+                                 weighted ? PERIODIC_PRIM_LS : PERIODIC_PRIM_ULS, *field.geometry, *field.config,
+                                 weighted, nodes->GetPrimitive(), 0, solver->GetnPrimVarGrad(), -1,
+                                 nodes->GetGradient_Primitive(), nodes->GetRmatrix());
+  }
+  END_SU2_OMP_PARALLEL
+  for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) {
+    CHECK(nodes->GetGradient_Primitive(i, 0, 0) == Approx(0).margin(1e-12));
+    CHECK(nodes->GetGradient_Primitive(i, 0, 1) == Approx(1).margin(1e-12));
+  }
+}
+
+TEST_CASE("Zone-local periodic pair intersections", "[Periodic][Geometry]") {
+  const bool absentPair = GENERATE(false, true);
+  UnitQuadTestCase field;
+  const auto start = field.config_options.find("MARKER_HEATFLUX=");
+  const auto end = field.config_options.find("VISCOSITY_MODEL=");
+  field.config_options.replace(
+      start, end - start,
+      absentPair ? "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n" : "MARKER_CUSTOM= (z_plus,z_minus)\n");
+  field.AddOption(absentPair ? "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0, "
+                               "other_minus,other_plus, 0,0,0, 0,0,0, 0,1,0)"
+                             : "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0, "
+                               "y_minus,y_plus, 0,0,0, 0,0,0, 0,1,0)");
+  field.InitConfig();
+  field.InitGeometry();
+  REQUIRE(field.config->GetnMarker_Periodic() == 4);
+  if (!absentPair && SU2_MPI::GetRank() != 0) {
+    /*--- Only rank zero has a local intersection; every rank must select the same communication path. ---*/
+    for (auto marker = 0u; marker < field.geometry->GetnMarker(); ++marker) {
+      if (field.config->GetMarker_All_PerBound(marker) != 2 && field.config->GetMarker_All_PerBound(marker) != 4)
+        continue;
+      for (auto v = 0ul; v < field.geometry->GetnVertex(marker); ++v)
+        field.geometry->nodes->SetVertex(field.geometry->vertex[marker][v]->GetNode(), -1, marker);
+    }
+  }
+  field.geometry->PreprocessPeriodicGeometry(*field.config);
+  CHECK(field.geometry->HasPeriodicPairIntersections() == !absentPair);
+}
+
+TEST_CASE("Invalid rotational intersection", "[.PeriodicInvalid]") {
+  UnitQuadTestCase field;
+  const auto start = field.config_options.find("MARKER_HEATFLUX=");
+  const auto end = field.config_options.find("VISCOSITY_MODEL=");
+  field.config_options.replace(start, end - start, "MARKER_CUSTOM= (z_plus,z_minus)\n");
+  field.AddOption(
+      "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,90, 1,0,0, "
+      "y_minus,y_plus, 0,0,0, 0,0,0, 0,1,0)");
+  field.InitConfig();
+  field.InitGeometry();
+  field.geometry->PreprocessPeriodicGeometry(*field.config);
+  FAIL("Rotational intersection was accepted");
+}
+
+TEST_CASE("Principal periodic rotation angle", "[Periodic][Geometry]") {
+  const su2double zero[3] = {};
+  CHECK(GeometryToolbox::RotationAngle(3, zero) == 0.0);
+  const su2double quarterTurn[3] = {0, 0, -PI_NUMBER / 2};
+  CHECK(GeometryToolbox::RotationAngle(2, quarterTurn) == Approx(PI_NUMBER / 2));
+  CHECK(GeometryToolbox::RotationAngle(3, quarterTurn) == Approx(PI_NUMBER / 2));
+  const su2double composed[3] = {PI_NUMBER / 2, PI_NUMBER / 2, 0};
+  CHECK(GeometryToolbox::RotationAngle(2, composed) == 0.0);
+  CHECK(GeometryToolbox::RotationAngle(3, composed) == Approx(2 * PI_NUMBER / 3));
+  const su2double halfTurn[3] = {PI_NUMBER, 0, 0};
+  CHECK(GeometryToolbox::RotationAngle(3, halfTurn) == Approx(PI_NUMBER));
+}
+
+TEST_CASE("Periodic edge weight cache rebuild", "[Periodic][Geometry]") {
+  UnitQuadTestCase field;
+  const auto start = field.config_options.find("MARKER_HEATFLUX=");
+  const auto end = field.config_options.find("VISCOSITY_MODEL=");
+  field.config_options.replace(start, end - start, "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n");
+  field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)");
+  field.InitConfig();
+  field.InitGeometry(true);
+  auto& geometry = *field.geometry;
+  geometry.PreprocessPeriodicGeometry(*field.config);
+  unsigned long edge = geometry.GetnEdge(), point = 0;
+  unsigned short marker = 0;
+  long vertex = -1;
+  for (auto iEdge = 0ul; iEdge < geometry.GetnEdge(); ++iEdge) {
+    if (geometry.GetPeriodicEdgeWeight(iEdge) == 1.0) continue;
+    const auto iPoint = geometry.edges->GetNode(iEdge, 0);
+    for (auto iMarker = 0u; iMarker < geometry.GetnMarker(); ++iMarker) {
+      if (field.config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) continue;
+      const auto iVertex = geometry.nodes->GetVertex(iPoint, iMarker);
+      if (iVertex < 0) continue;
+      edge = iEdge;
+      point = iPoint;
+      marker = iMarker;
+      vertex = iVertex;
+      break;
+    }
+    break;
+  }
+  const bool checked = edge < geometry.GetnEdge();
+  if (checked) {
+    CHECK(geometry.GetPeriodicEdgeWeight(edge) == 0.5);
+    geometry.nodes->SetVertex(point, -1, marker);
+  }
+  /*--- All ranks rebuild together, including those without a local candidate edge. ---*/
+  geometry.PreprocessPeriodicGeometry(*field.config);
+  if (checked) {
+    CHECK(geometry.GetPeriodicEdgeWeight(edge) == 1.0);
+    geometry.nodes->SetVertex(point, vertex, marker);
+  }
+  geometry.PreprocessPeriodicGeometry(*field.config);
+  if (checked) CHECK(geometry.GetPeriodicEdgeWeight(edge) == 0.5);
+  unsigned long local = checked, global = 0;
+  SU2_MPI::Allreduce(&local, &global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  REQUIRE(global > 0);
+}
+
+TEST_CASE("Coarse periodic boundary flags", "[Periodic]") {
+  UnitQuadTestCase field;
+  const auto start = field.config_options.find("MARKER_HEATFLUX=");
+  const auto end = field.config_options.find("VISCOSITY_MODEL=");
+  field.config_options.replace(start, end - start,
+                               "MARKER_HEATFLUX= (y_minus,0,y_plus,0)\nMARKER_CUSTOM= (z_plus,z_minus)\n");
+  field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)");
+  field.AddOption("MGLEVEL= 1");
+  field.InitConfig();
+  field.InitGeometry(true);
+  CMultiGridGeometry coarse(field.geometry.get(), field.config.get(), 1);
+  coarse.SetPoint_Connectivity(field.geometry.get());
+  coarse.SetVertex(field.geometry.get(), field.config.get());
+  unsigned long periodic = 0, missing = 0;
+  for (auto marker = 0u; marker < coarse.GetnMarker(); ++marker) {
+    if (field.config->GetMarker_All_KindBC(marker) == HEAT_FLUX) {
+      for (auto v = 0ul; v < coarse.GetnVertex(marker); ++v) {
+        const auto i = coarse.vertex[marker][v]->GetNode();
+        CHECK(coarse.nodes->GetPhysicalBoundary(i));
+        CHECK(coarse.nodes->GetSolidBoundary(i));
+        CHECK(coarse.nodes->GetViscousBoundary(i));
+      }
+    }
+    if (field.config->GetMarker_All_KindBC(marker) != PERIODIC_BOUNDARY) continue;
+    for (auto vertex = 0ul; vertex < coarse.GetnVertex(marker); ++vertex) {
+      ++periodic;
+      if (!coarse.nodes->GetPeriodicBoundary(coarse.vertex[marker][vertex]->GetNode())) ++missing;
+    }
+  }
+  unsigned long local[] = {periodic, missing}, global[2] = {};
+  SU2_MPI::Allreduce(local, global, 2, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  REQUIRE(global[0] > 0);
+  CHECK(global[1] == 0);
+}
+
+TEST_CASE("One-cell periodic stencil", "[Periodic][Gradients]") {
+  UnitQuadTestCase field;
+  const auto start = field.config_options.find("MARKER_HEATFLUX=");
+  const auto end = field.config_options.find("VISCOSITY_MODEL=");
+  field.config_options.replace(start, end - start, "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n");
+  field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)");
+  field.SetOption("MESH_BOX_SIZE= 2,5,5");
+  field.AddOption("NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES");
+  field.InitConfig();
+  field.InitGeometry(true);
+  field.geometry->MatchPeriodic(field.config.get(), 1);
+  field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get());
+  CHECK(field.geometry->HasPeriodicPairIntersections() == (field.config->GetnMarker_Periodic() > 2));
+  field.InitSolver();
+  auto* solver = field.solver[FLOW_SOL];
+  auto* nodes = solver->GetNodes();
+  for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) nodes->SetPrimitive(i, v, 3 + x[1] + x[2]);
+  }
+  solver->SetRotatePeriodic(false);
+  auto& gradient = nodes->GetGradient_Primitive();
+  auto& matrix = nodes->GetRmatrix();
+  SU2_OMP_PARALLEL {
+    computeGradientsLeastSquares(solver, MPI_QUANTITIES::PRIMITIVE_GRADIENT, PERIODIC_PRIM_LS, *field.geometry,
+                                 *field.config, true, nodes->GetPrimitive(), 0, solver->GetnPrimVarGrad(), -1, gradient,
+                                 matrix);
+  }
+  END_SU2_OMP_PARALLEL
+  unsigned long local = 0, global = 0;
+  for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    if (x[1] != 0.5 || x[2] != 0.5) continue;
+    ++local;
+    CHECK(field.geometry->nodes->GetnNeighbor(i) == 6);
+    CHECK(matrix(i, 0, 0) == Approx(2).margin(1e-12));
+    CHECK(gradient(i, 0, 1) == Approx(1).margin(1e-12));
+    CHECK(gradient(i, 0, 2) == Approx(1).margin(1e-12));
+  }
+  SU2_MPI::Allreduce(&local, &global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  REQUIRE(global == 2);
+}
diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp
index d8da75dc0960..d943c1696ee5 100644
--- a/UnitTests/UnitQuadTestCase.hpp
+++ b/UnitTests/UnitQuadTestCase.hpp
@@ -61,6 +61,16 @@ struct UnitQuadTestCase {
    */
   void AddOption(const std::string& optionLine) { config_options += optionLine + "\n"; }
 
+  /*! \brief Replace one existing base option without repeating its key. */
+  void SetOption(const std::string& optionLine) {
+    const auto key = optionLine.substr(0, optionLine.find('=') + 1);
+    const auto start = config_options.find(key);
+    if (start == std::string::npos)
+      AddOption(optionLine);
+    else
+      config_options.replace(start, config_options.find('\n', start) - start, optionLine);
+  }
+
   /*!
    * \brief Initialize the config structure
    */
@@ -83,10 +93,13 @@ struct UnitQuadTestCase {
   /*!
    * \brief Initialize the geometry
    */
-  void InitGeometry() {
+  void InitGeometry(bool partition = false) {
     cout.rdbuf(nullptr);
     {
-      auto aux_geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(config.get(), 0, 1));
+      const auto rank = partition ? SU2_MPI::GetRank() : 0;
+      const auto size = partition ? SU2_MPI::GetSize() : 1;
+      auto aux_geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(config.get(), rank, size));
+      if (partition) aux_geometry->SetColorGrid_Parallel(config.get());
       geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(aux_geometry.get(), config.get()));
     }
     geometry->SetSendReceive(config.get());
diff --git a/UnitTests/meson.build b/UnitTests/meson.build
index c6216606338f..8f67b4ef016f 100644
--- a/UnitTests/meson.build
+++ b/UnitTests/meson.build
@@ -53,6 +53,9 @@ if get_option('enable-tests')
         cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args]
     )
     test('Catch2 test driver', test_driver)
+    test('Reject intersecting rotational periodic pairs', python,
+         args: files('SU2_CFD/check_periodic_intersection_rejection.py') + [test_driver.full_path()],
+         timeout: 90)
     test('NEMO convergence field migration', python,
          args: files('SU2_CFD/output/check_nemo_convergence_migration.py') + [test_driver.full_path()],
          timeout: 300)