Fixes for rotational periodic boundaries: implicit solver, multigrid and limiters

Common/include/toolboxes/geometry_toolbox.hpp

@@ -27,6 +27,7 @@
 #pragma once
 
 #include <cmath>
+#include <algorithm>
 
 namespace GeometryToolbox {
 /// \addtogroup GeometryToolbox
@@ -215,6 +216,31 @@ inline void Rotate(const Scalar R[][nDim], const Scalar* O, const Scalar* d, Sca
   }
 }
 
+/*! \return Whether any of the three supplied rotation angles is nonzero. */
+template <class Scalar>
+inline bool HasRotation(const Scalar* angles) {
+  // Configured zero stays exact after degree conversion and negation; retain every nonzero rotation.
+  return angles[0] != 0.0 || angles[1] != 0.0 || angles[2] != 0.0;
+}
+
+/*! \brief Rotate component bounds in place, enclosing the rotated box. */
+template <class Scalar, int nDim>
+inline void RotateBox(const Scalar R[][nDim], Scalar* vMin, Scalar* vMax) {
+  using std::max;
+  using std::min;
+  Scalar rotMin[nDim] = {0.0}, rotMax[nDim] = {0.0};
+  for (int iDim = 0; iDim < nDim; ++iDim) {
+    for (int jDim = 0; jDim < nDim; ++jDim) {
+      const Scalar fromMin = R[iDim][jDim] * vMin[jDim];
+      const Scalar fromMax = R[iDim][jDim] * vMax[jDim];
+      rotMin[iDim] += min(fromMin, fromMax);
+      rotMax[iDim] += max(fromMin, fromMax);
+    }
+  }
+  std::copy_n(rotMin, nDim, vMin);
+  std::copy_n(rotMax, nDim, vMax);
+}
+
 /*! \brief Tangent projection  */
 template <class Mat, class Scalar, class Int>
 inline void TangentProjection(Int nDim, const Mat& tensor, const Scalar* vector, Scalar* proj) {

SU2_CFD/include/limiters/CLimiterDetails.hpp

@@ -72,6 +72,16 @@ struct LimiterHelpers
 {
   FORCEINLINE static Type epsilon() {return std::numeric_limits<passivedouble>::epsilon();}
 
+  /*! \brief MUSCL reconstruction increment using the displacement to the middle of an edge. */
+  template <class Int>
+  FORCEINLINE static Type reconstructionIncrement(Int nDim, const Type* halfEdge, const Type* gradient,
+                                                  const Type& value_i, const Type& value_j, const Type& kappa) {
+    Type proj = 0.0;
+    for (Int iDim = 0; iDim < nDim; ++iDim) proj += halfEdge[iDim] * gradient[iDim];
+    const Type cent = 0.5 * (value_j - value_i);
+    return umusclProjection(proj, cent, kappa);
+  }
+
   FORCEINLINE static Type umusclProjection(const Type& grad_proj, const Type& delta, const Type& kappa)
   {
     /*-------------------------------------------------------------------*/

SU2_CFD/include/limiters/computeLimiters_impl.hpp

@@ -109,11 +109,29 @@ void computeLimiters_impl(CSolver* solver,
 
   limiterDetails.preprocess(geometry, config, varBegin, varEnd, field);
 
+  /*--- With rotational periodicity, the first periodic comm. also brings the min/max
+   *    projections over the edges of the periodic matches (stored after each other),
+   *    because the limiters of the velocity cannot be compared across a rotation. ---*/
+
+  su2activematrix* periodicProj = nullptr;
+
+  if (periodic && (kindPeriodicComm1 == PERIODIC_LIM_PRIM_1))
+    periodicProj = solver->GetPeriodicProjections(geometry, config);
+
   /*--- Initialize all min/max field values if we have
    *    periodic comms. otherwise do it inside main loop. ---*/
 
   if (periodic)
   {
+    if (periodicProj != nullptr)
+    {
+      SU2_OMP_FOR_STAT(chunkSize)
+      for (auto iPoint = 0ul; iPoint < periodicProj->rows(); ++iPoint)
+        for (auto iVar = 0ul; iVar < periodicProj->cols(); ++iVar)
+          (*periodicProj)(iPoint,iVar) = 0.0;
+      END_SU2_OMP_FOR
+    }
+
     SU2_OMP_FOR_STAT(chunkSize)
     for (size_t iPoint = 0; iPoint < nPoint; ++iPoint)
       for (size_t iVar = varBegin; iVar < varEnd; ++iVar)
@@ -165,6 +183,23 @@ void computeLimiters_impl(CSolver* solver,
     for (size_t iVar = varBegin; iVar < varEnd; ++iVar)
       projMax[iVar] = projMin[iVar] = 0.0;
 
+    if (periodicProj != nullptr && nodes->GetPeriodicBoundary(iPoint))
+    {
+      /*--- Start from the min/max over the edges of the periodic matches. ---*/
+
+      const auto* projections = solver->GetPeriodicProjection(iPoint);
+      if (projections != nullptr) {
+        for (auto iVar = varBegin; iVar < varEnd; ++iVar) {
+          const auto& periodicMin = projections[iVar];
+          const auto& periodicMax = projections[periodicProj->cols()/2 + iVar];
+          AD::SetPreaccIn(periodicMin);
+          AD::SetPreaccIn(periodicMax);
+          projMin[iVar] = periodicMin;
+          projMax[iVar] = periodicMax;
+        }
+      }
+    }
+
     /*--- Compute max/min projection and values over direct neighbors. ---*/
 
     for (auto jPoint : geometry.nodes->GetPoints(iPoint)) {
@@ -175,22 +210,16 @@ void computeLimiters_impl(CSolver* solver,
       /*--- Distance vector from iPoint to face (middle of the edge). ---*/
 
       su2double dist_ij[nDim] = {0.0};
-
-      for(size_t iDim = 0; iDim < nDim; ++iDim)
+      for (size_t iDim = 0; iDim < nDim; ++iDim)
         dist_ij[iDim] = 0.5 * (coord_j[iDim] - coord_i[iDim]);
 
       /*--- Project each variable, update min/max. ---*/
 
       for(size_t iVar = varBegin; iVar < varEnd; ++iVar)
       {
-        su2double proj = 0.0;
-
-        for(size_t iDim = 0; iDim < nDim; ++iDim)
-          proj += dist_ij[iDim] * gradient(iPoint,iVar,iDim);
-
         AD::SetPreaccIn(field(jPoint,iVar));
-        const su2double cent = 0.5 * (field(jPoint,iVar) - field(iPoint,iVar));
-        proj = LimiterHelpers<>::umusclProjection(proj, cent, umusclKappa);
+        const su2double proj = LimiterHelpers<>::reconstructionIncrement(nDim, dist_ij,
+                              gradient[iPoint][iVar], field(iPoint,iVar), field(jPoint,iVar), umusclKappa);
 
         projMax[iVar] = max(projMax[iVar], proj);
         projMin[iVar] = min(projMin[iVar], proj);
@@ -224,7 +253,8 @@ void computeLimiters_impl(CSolver* solver,
   }
   END_SU2_OMP_FOR
 
-  /*--- Account for periodic effects, take the minimum limiter on each periodic pair. ---*/
+  /*--- Account for periodic effects, take the minimum limiter on each periodic pair
+   *    (except for the velocity with rotational periodicity). ---*/
   if (periodic)
   {
     for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic)

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -261,7 +261,7 @@ void CFVMFlowSolverBase<V, R>::CommunicateInitialState(CGeometry* geometry, cons
     CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
   }
   SetImplicitPeriodic(euler_implicit);
-  if (MGLevel == MESH_0) SetRotatePeriodic(true);
+  SetRotatePeriodic(true);
 
   /*--- Perform the MPI communication of the solution ---*/
 

SU2_CFD/include/solvers/CSolver.hpp

@@ -37,6 +37,7 @@
 #include <algorithm>
 #include <iostream>
 #include <set>
+#include <unordered_map>
 #include <stdlib.h>
 #include <stdio.h>
 
@@ -146,6 +147,8 @@ class CSolver {
 
   bool rotate_periodic;    /*!< \brief Flag that controls whether the periodic solution needs to be rotated for the solver. */
   bool implicit_periodic;  /*!< \brief Flag that controls whether the implicit system should be treated by the periodic BC comms. */
+  su2activematrix PeriodicProj;  /*!< \brief Min/max reconstruction increments at periodic receive points (for limiters). */
+  std::unordered_map<unsigned long, size_t> PeriodicProjIndex;  /*!< \brief Local point to compact projection row. */
 
   bool dynamic_grid;       /*!< \brief Flag that determines whether the grid is dynamic (moving or deforming + grid velocities). */
 
@@ -4231,6 +4234,22 @@ class CSolver {
    */
   inline void SetRotatePeriodic(bool val_rotate_periodic) { rotate_periodic = val_rotate_periodic; }
 
+  /*!
+   * \brief Storage for the limiters with rotational periodicity: the min and max, over the edges of the periodic
+   *        matches of each point, of the reconstruction increments (communicated with PERIODIC_LIM_PRIM_1).
+   * \param[in] geometry - Periodic receive points of this mesh level.
+   * \param[in] config - Definition of the particular problem.
+   * \return The matrix (unique periodic receive points x 2*nPrimVarGrad, min then max), nullptr without rotation.
+   */
+  su2activematrix* GetPeriodicProjections(const CGeometry& geometry, const CConfig& config);
+
+  /*! \brief Reconstruction increment bounds for a periodic receive point, nullptr for other points. */
+  inline su2double* GetPeriodicProjection(unsigned long iPoint) {
+    const auto& indices = PeriodicProjIndex;
+    const auto it = indices.find(iPoint);
+    return it == indices.end() ? nullptr : PeriodicProj[it->second];
+  }
+
   /*!
    * \brief Retrieve the solver name for output purposes.
    * \returns Name of the solver.