Fix periodic neighbor counts, intersecting stencils and coarse flags

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.

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>

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++) {

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));
       }
     }
 

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 ---*/
 

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);
         }
       }
 

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);
 

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. */
 

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

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  ---*/