Fix periodic neighbor counts, intersecting stencils and coarse flags

Common/include/geometry/CGeometry.hpp

@@ -62,6 +62,7 @@ extern "C" {
 
 #include "../CConfig.hpp"
 #include "../toolboxes/graph_toolbox.hpp"
+#include "../toolboxes/geometry_toolbox.hpp"
 #include "../adt/CADTElemClass.hpp"
 
 using namespace std;
@@ -517,6 +518,32 @@ class CGeometry {
    */
   void AllocatePeriodicComms(unsigned short val_countPerPeriodicPoint);
 
+  /*! \brief Fraction of an edge stencil owned by this periodic copy. */
+  inline passivedouble GetPeriodicEdgeWeight(unsigned long iPoint, unsigned long jPoint, const CConfig& config) const {
+    if (!nodes->GetPeriodicBoundary(iPoint) || !nodes->GetPeriodicBoundary(jPoint)) return 1.0;
+    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);
+        su2double rotation[3][3];
+        GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation);
+        const auto cosAngle = SU2_TYPE::GetValue(0.5 * (rotation[0][0] + rotation[1][1] + rotation[2][2] - 1));
+        weight *= acos(std::max(-1.0, std::min(1.0, cosAngle))) / (2 * PI_NUMBER);
+      } else {
+        weight *= 0.5;
+      }
+    }
+    return weight;
+  }
+
   /*!
    * \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/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

@@ -253,6 +253,8 @@ void computeGradientsLeastSquares(CSolver* solver,
       if (weight > 0.0)
       {
         weight = 1.0 / weight;
+        if (config.GetnMarker_Periodic() > 2)
+          weight *= geometry.GetPeriodicEdgeWeight(iPoint, jPoint, config);
 
         for (size_t iDim = 0; iDim < nDim; ++iDim)
           for (size_t jDim = iDim; jDim < nDim; ++jDim)

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -807,8 +807,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 = config->GetnMarker_Periodic() > 2 ?
+                                  geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) : 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);
   if (MGLevel == MESH_0) SetRotatePeriodic(true);
 

SU2_CFD/include/solvers/CSolver.hpp

@@ -133,6 +133,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

@@ -2404,13 +2404,16 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *co
       /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/
       if (boundary_i && !boundary_j) continue;
 
+      const auto weight = config->GetnMarker_Periodic() > 2 ?
+                              geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) : 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  ---*/

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;
@@ -352,6 +352,18 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
     SU2_MPI::Error("The NEMO solvers do not support rotational periodicity yet.", CURRENT_FUNCTION);
   }
 
+  const bool intersectingPairs = config->GetnMarker_Periodic() > 2;
+  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 jPoint : geometry->nodes->GetPoints(iPoint))
+        periodicNeighborCount(iPoint) += geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config);
+    }
+    END_SU2_OMP_FOR
+  }
+
   /*--- Local variables ---*/
 
   bool boundary_i, boundary_j;
@@ -413,7 +425,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. ---*/
 
@@ -501,20 +512,25 @@ 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)) {
 
               /*--- 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(iPoint, jPoint, *config) == 1.0)
                 nNeighbor++;
             }
 
             /*--- Store the number of neighbors in bufffer. ---*/
 
-            bufSSend[buf_offset] = nNeighbor;
+            bufDSend[buf_offset] = nNeighbor;
 
             break;
 
@@ -605,6 +621,16 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
           case PERIODIC_LAPLACIAN:
 
+            if (intersectingPairs) {
+              for (auto iField = 0u; iField < nVar; ++iField)
+                bufDSend[buf_offset + iField] = base_nodes->GetUndivided_Laplacian(iPoint, iField);
+              if (rotate_periodic) Rotate(zeros, &bufDSend[buf_offset + 1], &Und_Lapl[1]);
+              if (rotate_periodic)
+                for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate)
+                  bufDSend[buf_offset + 1 + iCoordinate] = Und_Lapl[1 + iCoordinate];
+              break;
+            }
+
             /*--- For JST, the undivided Laplacian must be computed
              consistently by using the complete control volume info
              from both sides of the periodic face. ---*/
@@ -617,7 +643,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
               /*--- 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(iPoint, jPoint, *config) == 1.0) {
 
                 /*--- Solution differences ---*/
 
@@ -671,6 +697,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
@@ -683,7 +714,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
               /*--- 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(iPoint, jPoint, *config) < 1.0) continue;
 
               /*--- Use density instead of pressure for incomp. flows. ---*/
 
@@ -703,7 +734,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);
@@ -781,6 +812,42 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
             /*--- Set a flag for unweighted or weighted least-squares. ---*/
 
+            if (intersectingPairs) {
+              /*--- Rmatrix stores the upper triangle of the normal equations;
+               * entry (2,1) is a duplicate of (0,2) used by the LS factorization. ---*/
+              su2double matrix[3][3] = {}, rotated[3][3] = {};
+              for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate)
+                for (auto jDim = 0u; jDim < nDim; ++jDim)
+                  matrix[iCoordinate][jDim] =
+                      base_nodes->GetRmatrix(iPoint, min(iCoordinate,jDim), max(iCoordinate,jDim));
+              for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate)
+                for (auto jDim = 0u; jDim < nDim; ++jDim)
+                  for (auto kDim = 0u; kDim < nDim; ++kDim)
+                    for (auto lDim = 0u; lDim < nDim; ++lDim) {
+                      const auto qik = nDim == 2 ? rotMatrix2D[iCoordinate][kDim] : rotMatrix3D[iCoordinate][kDim];
+                      const auto qjl = nDim == 2 ? rotMatrix2D[jDim][lDim] : rotMatrix3D[jDim][lDim];
+                      rotated[iCoordinate][jDim] += qik * matrix[kDim][lDim] * qjl;
+                    }
+              for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate)
+                for (auto jDim = 0u; jDim < nDim; ++jDim)
+                  bufDSend[buf_offset++] = iCoordinate <= jDim ? rotated[iCoordinate][jDim] :
+                      (nDim == 3 && iCoordinate == 2 && jDim == 1 ? rotated[0][2] : su2double(0));
+              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:
@@ -826,7 +893,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
               /*--- 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(iPoint, jPoint, *config) == 1.0) {
 
                 /*--- Get coordinates for the neighbor point. ---*/
 
@@ -1041,7 +1108,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. ---*/
 
@@ -1118,13 +1184,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 (config->GetnMarker_Periodic() > 2) {
+                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:
@@ -2280,6 +2346,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 (config->GetnMarker_Periodic() > 2) delta *= geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config);
         base_nodes->AddUnd_Lapl(iPoint, iVar, delta);
       }
     }