[WIP] Search translated wall images when computing periodic wall distance

Common/include/adt/CADTElemClass.hpp

@@ -31,6 +31,8 @@
 #include "./CBBoxTargetClass.hpp"
 #include "../parallelization/omp_structure.hpp"
 
+class CConfig;
+
 /*!
  * \class CADTElemClass
  * \ingroup ADT
@@ -42,6 +44,9 @@ class CADTElemClass : public CADTBaseClass {
  private:
   unsigned short nDim; /*!< \brief Number of spatial dimensions. */
 
+  vector<array<su2double, 3>> periodicTranslations, periodicDual;
+  array<su2double, 3> wallMin{}, wallMax{};
+
   vector<su2double> coorPoints; /*!< \brief Vector, which contains the coordinates
                                             of the points in the ADT. */
   vector<su2double> BBoxCoor;   /*!< \brief Vector, which contains the coordinates
@@ -61,8 +66,8 @@ class CADTElemClass : public CADTBaseClass {
   vector<int> ranksOfElems;            /*!< \brief Vector, which contains the ranks
                                                    of the elements in the ADT. */
 #ifdef HAVE_OMP
-  vector<vector<CBBoxTargetClass> > BBoxTargets; /*!< \brief Vector, used to store possible bounding box
-                                                             candidates during the nearest element search. */
+  vector<vector<CBBoxTargetClass>> BBoxTargets; /*!< \brief Vector, used to store possible bounding box
+                                                            candidates during the nearest element search. */
 #else
   array<vector<CBBoxTargetClass>, 1> BBoxTargets;
 #endif
@@ -84,6 +89,9 @@ class CADTElemClass : public CADTBaseClass {
                 vector<unsigned short>& val_VTKElem, vector<unsigned short>& val_markerID,
                 vector<unsigned long>& val_elemID, const bool globalTree);
 
+  /*! \brief Enable exact translated wall-image searches for this wall source zone. */
+  void SetPeriodicWallSearch(const CConfig* config);
+
   /*!
    * \brief Function, which determines the element that contains the given coordinate.
    * \note This simply forwards the call to the implementation function selecting the right
@@ -120,9 +128,13 @@ class CADTElemClass : public CADTBaseClass {
     const auto iThread = omp_get_thread_num();
     DetermineNearestElement_impl(BBoxTargets[iThread], FrontLeaves[iThread], FrontLeavesNew[iThread], coor, dist,
                                  markerID, elemID, rankID);
+    if (!periodicTranslations.empty()) DetermineNearestPeriodicElement(coor, dist, markerID, elemID, rankID);
   }
 
  private:
+  void DetermineNearestPeriodicElement(const su2double* coor, su2double& dist, unsigned short& markerID,
+                                       unsigned long& elemID, int& rankID);
+
   /*!
    * \brief Implementation of DetermineContainingElement.
    * \note Working variables (first two) passed explicitly for thread safety.

Common/src/adt/CADTElemClass.cpp

@@ -26,6 +26,9 @@
  */
 
 #include "../../include/adt/CADTElemClass.hpp"
+#include "../../include/CConfig.hpp"
+#include "../../include/linear_algebra/blas_structure.hpp"
+#include "../../include/toolboxes/geometry_toolbox.hpp"
 #include "../../include/parallelization/mpi_structure.hpp"
 #include "../../include/option_structure.hpp"
 
@@ -260,6 +263,121 @@ CADTElemClass::CADTElemClass(unsigned short val_nDim, vector<su2double>& val_coo
   for (auto& vec : FrontLeavesNew) vec.reserve(200);
 }
 
+void CADTElemClass::SetPeriodicWallSearch(const CConfig* config) {
+  periodicTranslations.clear();
+  periodicDual.clear();
+  if (IsEmpty() || !config->GetnMarker_Periodic()) return;
+  for (auto marker = 0u; marker < config->GetnMarker_CfgFile(); ++marker) {
+    const auto tag = config->GetMarker_CfgFile_TagBound(marker);
+    if (config->GetMarker_CfgFile_KindBC(tag) != PERIODIC_BOUNDARY ||
+        config->GetMarker_CfgFile_PerBound(tag) > config->GetnMarker_Periodic() / 2)
+      continue;
+    const auto* angles = config->GetPeriodicRotAngles(tag);
+    su2double rotation[3][3];
+    GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation);
+    for (auto i = 0u; i < nDim; ++i)
+      for (auto j = 0u; j < nDim; ++j)
+        if (fabs(rotation[i][j] - (i == j)) > 1e-12) {
+          periodicTranslations.clear();
+          return;
+        }
+    array<su2double, 3> translation{};
+    const auto* shift = config->GetPeriodicTranslation(tag);
+    std::copy(shift, shift + nDim, translation.begin());
+    periodicTranslations.push_back(translation);
+  }
+  const auto count = periodicTranslations.size();
+  if (!count || count > nDim) {
+    periodicTranslations.clear();
+    return;
+  }
+  // ponytail: independent translations; dependent generators need lattice reduction.
+  vector<array<su2double, 3>> orthogonal;
+  for (const auto& translation : periodicTranslations) {
+    auto vector = translation;
+    for (const auto& basis : orthogonal) {
+      const auto projection = GeometryToolbox::DotProduct(nDim, vector.data(), basis.data());
+      for (auto d = 0u; d < nDim; ++d) vector[d] -= projection * basis[d];
+    }
+    const auto length = GeometryToolbox::Norm(nDim, vector.data());
+    if (length <= 1e-12 * GeometryToolbox::Norm(nDim, translation.data())) {
+      periodicTranslations.clear();
+      return;
+    }
+    for (auto d = 0u; d < nDim; ++d) vector[d] /= length;
+    orthogonal.push_back(vector);
+  }
+  su2activematrix gram(count, count);
+  for (auto i = 0u; i < count; ++i)
+    for (auto j = 0u; j < count; ++j)
+      gram(i, j) = GeometryToolbox::DotProduct(nDim, periodicTranslations[i].data(), periodicTranslations[j].data());
+  /*--- Distinct translation pairs form the basis of the periodic cell. ---*/
+  CBlasStructure::inverse(count, gram);
+  periodicDual.resize(count);
+  for (auto i = 0u; i < count; ++i)
+    for (auto d = 0u; d < nDim; ++d)
+      for (auto j = 0u; j < count; ++j) periodicDual[i][d] += gram(i, j) * periodicTranslations[j][d];
+  for (auto d = 0u; d < nDim; ++d) {
+    wallMin[d] = wallMax[d] = coorPoints[d];
+    for (auto i = d; i < coorPoints.size(); i += nDim) {
+      wallMin[d] = min(wallMin[d], coorPoints[i]);
+      wallMax[d] = max(wallMax[d], coorPoints[i]);
+    }
+  }
+}
+
+void CADTElemClass::DetermineNearestPeriodicElement(const su2double* coor, su2double& dist, unsigned short& markerID,
+                                                    unsigned long& elemID, int& rankID) {
+  /*--- An improving image must lie inside the wall bounding box expanded by
+   * the current distance. Project that box onto the dual lattice basis to
+   * bound every integer shift, including skew cells and images beyond +/-1. ---*/
+  array<long, 3> lower{}, upper{}, index{};
+  const auto count = periodicTranslations.size();
+  for (auto i = 0u; i < count; ++i) {
+    su2double lo = 0, hi = 0;
+    for (auto d = 0u; d < nDim; ++d) {
+      const auto a = periodicDual[i][d] * (wallMin[d] - coor[d]);
+      const auto b = periodicDual[i][d] * (wallMax[d] - coor[d]);
+      lo += min(a, b);
+      hi += max(a, b);
+    }
+    const auto radius = dist * GeometryToolbox::Norm(nDim, periodicDual[i].data());
+    lower[i] = static_cast<long>(ceil(SU2_TYPE::GetValue(lo - radius)));
+    upper[i] = static_cast<long>(floor(SU2_TYPE::GetValue(hi + radius)));
+    if (lower[i] > upper[i]) return;
+    index[i] = lower[i];
+  }
+  const auto thread = omp_get_thread_num();
+  for (;;) {
+    array<su2double, 3> query{};
+    for (auto d = 0u; d < nDim; ++d) {
+      query[d] = coor[d];
+      for (auto i = 0u; i < count; ++i) query[d] += index[i] * periodicTranslations[i][d];
+    }
+    su2double candidate;
+    unsigned short marker;
+    unsigned long elem;
+    int rank;
+    DetermineNearestElement_impl(BBoxTargets[thread], FrontLeaves[thread], FrontLeavesNew[thread], query.data(),
+                                 candidate, marker, elem, rank);
+    if (candidate < dist) {
+      dist = candidate;
+      markerID = marker;
+      elemID = elem;
+      rankID = rank;
+    }
+    auto i = 0u;
+    for (; i < count; ++i) {
+      if (index[i] < upper[i]) {
+        ++index[i];
+        break;
+      }
+      index[i] = lower[i];
+    }
+    if (i == count) break;
+  }
+}
+
 bool CADTElemClass::DetermineContainingElement_impl(vector<unsigned long>& frontLeaves,
                                                     vector<unsigned long>& frontLeavesNew, const su2double* coor,
                                                     unsigned short& markerID, unsigned long& elemID, int& rankID,

Common/src/fem/fem_wall_distance.cpp

@@ -113,6 +113,7 @@ std::unique_ptr<CADTElemClass> CMeshFEM_DG::ComputeViscousWallADT(const CConfig*
   std::unique_ptr<CADTElemClass> WallADT(
       new CADTElemClass(nDim, surfaceCoor, surfaceConn, VTK_TypeElem, markerIDs, elemIDs, true));
 
+  WallADT->SetPeriodicWallSearch(config);
   return WallADT;
 }
 

Common/src/geometry/CPhysicalGeometry.cpp

@@ -10252,6 +10252,7 @@ std::unique_ptr<CADTElemClass> CPhysicalGeometry::ComputeViscousWallADT(const CC
   std::unique_ptr<CADTElemClass> WallADT(
       new CADTElemClass(nDim, surfaceCoor, surfaceConn, VTK_TypeElem, markerIDs, elemIDs, true));
 
+  WallADT->SetPeriodicWallSearch(config);
   return WallADT;
 }
 

UnitTests/Common/geometry/CGeometry_test.cpp

@@ -140,3 +140,49 @@ TEST_CASE("Set bound control volume", "[Geometry]") {
   CHECK(TestCase->geometry->vertex[3][2]->GetNormal()[1] == -0.0625);
   CHECK(TestCase->geometry->vertex[5][3]->GetNormal()[2] == 0.03125);
 }
+
+TEST_CASE("Periodic wall distance across translated images", "[Periodic][WallDistance]") {
+  UnitQuadTestCase field;
+  field.SetOption("SOLVER= RANS");
+  field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION");
+  field.SetOption("MARKER_CUSTOM= (z_plus,z_minus)");
+  field.AddOption("KIND_TURB_MODEL= SA");
+  field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)");
+  field.InitConfig();
+  field.InitGeometry(true);
+  for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    const auto hill = fabs(x[0] - 0.75) < 1e-12 ? 0.4 : 0.0;
+    field.geometry->nodes->SetCoord(i, 1, hill + (1 - hill) * x[1]);
+  }
+  field.geometry->SetWallDistance(std::numeric_limits<su2double>::max());
+  auto wall = field.geometry->ComputeViscousWallADT(field.config.get());
+  field.geometry->SetWallDistance(wall.get(), field.config.get(), 0);
+  su2double left = 0, right = 0, image = 0;
+  unsigned long checked = 0;
+  for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    if (fabs(x[1] - 0.5) > 1e-12 || fabs(x[2] - 0.5) > 1e-12) continue;
+    if (fabs(x[0]) < 1e-12) {
+      left = field.geometry->nodes->GetWall_Distance(i);
+      su2double shifted[] = {x[0] + 1, x[1], x[2]};
+      unsigned short marker;
+      unsigned long elem;
+      int rank;
+      wall->DetermineNearestElement(shifted, image, marker, elem, rank);
+      ++checked;
+    }
+    if (fabs(x[0] - 1) < 1e-12) right = field.geometry->nodes->GetWall_Distance(i);
+  }
+  su2double globalLeft = 0, globalRight = 0, globalImage = 0;
+  unsigned long total = 0;
+  SU2_MPI::Allreduce(&left, &globalLeft, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&right, &globalRight, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&image, &globalImage, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&checked, &total, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  REQUIRE(total == 1);
+  const auto expected = 0.5 / sqrt(1.0 + 1.6 * 1.6);
+  CHECK(globalImage == Approx(expected));
+  CHECK(globalRight == Approx(expected));
+  CHECK(globalLeft == Approx(expected));
+}

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