From 811dcb00325c87bb2e25bfe007487d15b7b272c2 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:44:30 +0200 Subject: [PATCH] Search independent translated wall images in periodic wall-distance queries --- Common/include/adt/CADTElemClass.hpp | 16 ++- Common/src/adt/CADTElemClass.cpp | 118 +++++++++++++++++++ Common/src/fem/fem_wall_distance.cpp | 1 + Common/src/geometry/CPhysicalGeometry.cpp | 1 + UnitTests/Common/geometry/CGeometry_test.cpp | 46 ++++++++ UnitTests/UnitQuadTestCase.hpp | 17 ++- 6 files changed, 195 insertions(+), 4 deletions(-) diff --git a/Common/include/adt/CADTElemClass.hpp b/Common/include/adt/CADTElemClass.hpp index 9222a1e0fd80..9e0986191c4f 100644 --- a/Common/include/adt/CADTElemClass.hpp +++ b/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> periodicTranslations, periodicDual; + array wallMin{}, wallMax{}; + vector coorPoints; /*!< \brief Vector, which contains the coordinates of the points in the ADT. */ vector BBoxCoor; /*!< \brief Vector, which contains the coordinates @@ -61,8 +66,8 @@ class CADTElemClass : public CADTBaseClass { vector ranksOfElems; /*!< \brief Vector, which contains the ranks of the elements in the ADT. */ #ifdef HAVE_OMP - vector > BBoxTargets; /*!< \brief Vector, used to store possible bounding box - candidates during the nearest element search. */ + vector> BBoxTargets; /*!< \brief Vector, used to store possible bounding box + candidates during the nearest element search. */ #else array, 1> BBoxTargets; #endif @@ -84,6 +89,9 @@ class CADTElemClass : public CADTBaseClass { vector& val_VTKElem, vector& val_markerID, vector& 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. diff --git a/Common/src/adt/CADTElemClass.cpp b/Common/src/adt/CADTElemClass.cpp index f65d9115edec..e49ae9fe1590 100644 --- a/Common/src/adt/CADTElemClass.cpp +++ b/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& 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 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> 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 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(ceil(SU2_TYPE::GetValue(lo - radius))); + upper[i] = static_cast(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 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& frontLeaves, vector& frontLeavesNew, const su2double* coor, unsigned short& markerID, unsigned long& elemID, int& rankID, diff --git a/Common/src/fem/fem_wall_distance.cpp b/Common/src/fem/fem_wall_distance.cpp index 022707c6868a..7c63090f29a2 100644 --- a/Common/src/fem/fem_wall_distance.cpp +++ b/Common/src/fem/fem_wall_distance.cpp @@ -113,6 +113,7 @@ std::unique_ptr CMeshFEM_DG::ComputeViscousWallADT(const CConfig* std::unique_ptr WallADT( new CADTElemClass(nDim, surfaceCoor, surfaceConn, VTK_TypeElem, markerIDs, elemIDs, true)); + WallADT->SetPeriodicWallSearch(config); return WallADT; } diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 6928c2d922b5..1f117ca34442 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -10252,6 +10252,7 @@ std::unique_ptr CPhysicalGeometry::ComputeViscousWallADT(const CC std::unique_ptr WallADT( new CADTElemClass(nDim, surfaceCoor, surfaceConn, VTK_TypeElem, markerIDs, elemIDs, true)); + WallADT->SetPeriodicWallSearch(config); return WallADT; } diff --git a/UnitTests/Common/geometry/CGeometry_test.cpp b/UnitTests/Common/geometry/CGeometry_test.cpp index 22556b16055d..f004dcb00a60 100644 --- a/UnitTests/Common/geometry/CGeometry_test.cpp +++ b/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::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)); +} 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(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(new CPhysicalGeometry(config.get(), rank, size)); + if (partition) aux_geometry->SetColorGrid_Parallel(config.get()); geometry = std::unique_ptr(new CPhysicalGeometry(aux_geometry.get(), config.get())); } geometry->SetSendReceive(config.get());