From 2990df482dee563cbb910294b9e99d752b01cc9b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:32 +0200 Subject: [PATCH 01/11] Rotate the columns of the Jacobian block at rotational periodic boundaries The residual of the periodic match is added as Q*R, and the solution of the match is Q*U, so the block added to the diagonal must be Q*J*Q^T. Only the rows were rotated. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- SU2_CFD/src/solvers/CSolver.cpp | 11 ++- .../navierstokes/periodic2D/no_limiter.cfg | 73 +++++++++++++++++++ TestCases/parallel_regression.py | 14 +++- TestCases/serial_regression.py | 6 +- 4 files changed, 97 insertions(+), 7 deletions(-) create mode 100644 TestCases/navierstokes/periodic2D/no_limiter.cfg diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 56cee4b6ad52..2f62f3241f36 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -553,7 +553,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, } } - /*--- Rotate the momentum columns of the Jacobian. ---*/ + /*--- Rotate the momentum rows and columns of the Jacobian, the residual of the + periodic match is Q*R(Q^T*U), so its Jacobian is Q*J*Q^T. First the rows. ---*/ if (rotate_periodic) { for (iVar = 0; iVar < nVar; iVar++) { @@ -569,6 +570,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, rotMatrix3D[2][2]*block(3, iVar); } } + + /*--- Then the columns, i.e. the momentum part of each row. ---*/ + + for (auto iRow = 0u; iRow < nVar; iRow++) { + su2double rotated[3] = {0.0}; + Rotate(zeros, &jacBlock[iRow][1], rotated); + for (auto jDim = 0u; jDim < nDim; jDim++) jacBlock[iRow][1+jDim] = rotated[jDim]; + } } /*--- Load the Jacobian terms into the buffer for sending. ---*/ diff --git a/TestCases/navierstokes/periodic2D/no_limiter.cfg b/TestCases/navierstokes/periodic2D/no_limiter.cfg new file mode 100644 index 000000000000..1c0a9320b94b --- /dev/null +++ b/TestCases/navierstokes/periodic2D/no_limiter.cfg @@ -0,0 +1,73 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Same as config.cfg without slope limiter (tests the Jacobian of the periodic points). +% +SOLVER= NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 0.1 +AOA= 22.5 +INIT_OPTION= TD_CONDITIONS +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_TEMPERATURE= 300 +FREESTREAM_PRESSURE= 99000 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.00 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1 +REF_AREA= 1 +% +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.87 +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.001 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_PERIODIC= ( per1, per2, 0,0,0, 0,0,45, 0,0,0 ) +MARKER_OUTLET= ( outlet, 99000.0 ) +INLET_TYPE= TOTAL_CONDITIONS +MARKER_INLET= ( inlet, 300, 100000, 1, 0, 0 ) +SPECIFIED_INLET_PROFILE= YES +INLET_FILENAME= inlet.dat +% +MARKER_ANALYZE= ( outlet, inlet ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 100 +CFL_ADAPT= NO +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 0.1 +LINEAR_SOLVER_ITER= 4 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -11 +CONV_STARTITER= 0 +ITER= 5000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= sector.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_RES, LINSOL_RESIDUAL, SURFACE_PRESSURE_DROP ) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index e85b5b7be660..8ef9648ae23b 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -421,6 +421,14 @@ def main(): poiseuille_profile.tol = [0.001, 0.001, 1e-5, 1e-5, 1e-5] test_list.append(poiseuille_profile) + # 2D rotational periodic sector without limiter (Jacobian of the periodic points) + periodic2d_no_limiter = TestCase('periodic2d_no_limiter') + periodic2d_no_limiter.cfg_dir = "navierstokes/periodic2D" + periodic2d_no_limiter.cfg_file = "no_limiter.cfg" + periodic2d_no_limiter.test_iter = 100 + periodic2d_no_limiter.test_vals = [-3.216881, -0.582813, -0.629049, 2.266141, -1.056877, -811.470000] + test_list.append(periodic2d_no_limiter) + ########################## ### Compressible RANS ### ########################## @@ -1244,7 +1252,7 @@ def main(): Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 Aachen_3D_restart.tol = 0.00001 - Aachen_3D_restart.test_vals = [-7.701420, -8.504728, -6.014939, -6.468223, -5.801124, -4.607179, -5.550665, -5.300778, -3.804188, -5.255983, -5.763060, -3.609605, -2.229249, -2.880453, -0.563469] + Aachen_3D_restart.test_vals = [-7.701423, -8.504851, -6.014951, -6.472219, -5.802012, -4.609768, -5.550659, -5.300718, -3.804222, -5.255983, -5.763064, -3.609604, -2.229253, -2.880518, -0.563484] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart @@ -1252,7 +1260,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.941917, -12.212515, -19.254664, -13.545311, -19.087161, -13.454459, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.910586, -12.203941, -19.198816, -13.489664, -19.033845, -13.400589, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -1277,7 +1285,7 @@ def main(): multi_interface.cfg_dir = "turbomachinery/multi_interface" multi_interface.cfg_file = "multi_interface_rst.cfg" multi_interface.test_iter = 5 - multi_interface.test_vals = [-8.632227, -8.894736, -9.348706] + multi_interface.test_vals = [-8.634558, -8.895554, -9.348754] multi_interface.test_vals_aarch64 = [-8.632227, -8.894736, -9.348706] test_list.append(multi_interface) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index c34fe6f47225..e1a8c131e34e 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -949,7 +949,7 @@ def main(): Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 - Aachen_3D_restart.test_vals = [-7.701421, -8.504727, -6.014939, -6.468221, -5.801125, -4.607173, -5.550665, -5.300779, -3.804187, -5.255982, -5.763060, -3.609601, -2.229250, -2.880453, -0.563470] + Aachen_3D_restart.test_vals = [-7.701424, -8.504850, -6.014951, -6.472217, -5.802013, -4.609762, -5.550659, -5.300718, -3.804222, -5.255982, -5.763064, -3.609600, -2.229254, -2.880518, -0.563485] Aachen_3D_restart.enabled_with_asan = False test_list.append(Aachen_3D_restart) @@ -958,7 +958,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.944235, -12.212620, -19.261137, -13.549357, -19.083828, -13.444697, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.916186, -12.202434, -19.225637, -13.511423, -19.026840, -13.394199, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -983,7 +983,7 @@ def main(): multi_interface.cfg_dir = "turbomachinery/multi_interface" multi_interface.cfg_file = "multi_interface_rst.cfg" multi_interface.test_iter = 5 - multi_interface.test_vals = [-8.632227, -8.894736, -9.348706] + multi_interface.test_vals = [-8.634558, -8.895554, -9.348754] multi_interface.test_vals_aarch64 = [-8.632227, -8.894736, -9.348706] test_list.append(multi_interface) From ce585c742212d70f453bb358578cf8b525a92b84 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:32 +0200 Subject: [PATCH 02/11] Rotate the periodic solution on all multigrid levels The coarse levels summed the momentum residuals, gradients and Jacobian blocks of a rotational periodic pair without rotating them. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- .../include/solvers/CFVMFlowSolverBase.inl | 2 +- .../navierstokes/periodic2D/multigrid.cfg | 78 +++++++++++++++++++ TestCases/parallel_regression.py | 8 ++ 3 files changed, 87 insertions(+), 1 deletion(-) create mode 100644 TestCases/navierstokes/periodic2D/multigrid.cfg diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index cb44757c6280..3a46fc8c4db6 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -261,7 +261,7 @@ void CFVMFlowSolverBase::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 ---*/ diff --git a/TestCases/navierstokes/periodic2D/multigrid.cfg b/TestCases/navierstokes/periodic2D/multigrid.cfg new file mode 100644 index 000000000000..72c8e9933439 --- /dev/null +++ b/TestCases/navierstokes/periodic2D/multigrid.cfg @@ -0,0 +1,78 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Same as config.cfg with multigrid, without slope limiter, lower CFL. +% +SOLVER= NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 0.1 +AOA= 22.5 +INIT_OPTION= TD_CONDITIONS +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_TEMPERATURE= 300 +FREESTREAM_PRESSURE= 99000 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.00 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1 +REF_AREA= 1 +% +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.87 +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.001 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_PERIODIC= ( per1, per2, 0,0,0, 0,0,45, 0,0,0 ) +MARKER_OUTLET= ( outlet, 99000.0 ) +INLET_TYPE= TOTAL_CONDITIONS +MARKER_INLET= ( inlet, 300, 100000, 1, 0, 0 ) +SPECIFIED_INLET_PROFILE= YES +INLET_FILENAME= inlet.dat +% +MARKER_ANALYZE= ( outlet, inlet ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 20 +CFL_ADAPT= NO +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 0.1 +LINEAR_SOLVER_ITER= 4 + +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 2 +MG_MIN_MESHSIZE= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -11 +CONV_STARTITER= 0 +ITER= 5000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= sector.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_RES, LINSOL_RESIDUAL, SURFACE_PRESSURE_DROP ) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 8ef9648ae23b..5bf6f7fb81f3 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -429,6 +429,14 @@ def main(): periodic2d_no_limiter.test_vals = [-3.216881, -0.582813, -0.629049, 2.266141, -1.056877, -811.470000] test_list.append(periodic2d_no_limiter) + # 2D rotational periodic sector with multigrid + periodic2d_multigrid = TestCase('periodic2d_multigrid') + periodic2d_multigrid.cfg_dir = "navierstokes/periodic2D" + periodic2d_multigrid.cfg_file = "multigrid.cfg" + periodic2d_multigrid.test_iter = 300 + periodic2d_multigrid.test_vals = [-4.422778, -1.615251, -1.453736, 1.047350, -1.410225, -2060.600000] + test_list.append(periodic2d_multigrid) + ########################## ### Compressible RANS ### ########################## From c39428c1daf24b2ab4d7f7ab17b9dc1e8cd9a083 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:32 +0200 Subject: [PATCH 03/11] Compute the limiters of rotational periodic points with the complete stencil The limiters of the velocity components were rotated like a vector when taking the minimum over a periodic pair, and the min and max velocity vectors were rotated instead of the velocity of each neighbour. Now each side sends, in the frame of its match, the min/max of the rotated velocities and of the rotated reconstruction increments, and the limiter is computed once. Nothing changes without rotation. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- .../include/limiters/computeLimiters_impl.hpp | 36 ++++- SU2_CFD/include/solvers/CSolver.hpp | 9 ++ SU2_CFD/src/solvers/CSolver.cpp | 143 ++++++++++++++++-- TestCases/hybrid_regression.py | 6 +- TestCases/parallel_regression.py | 20 ++- 5 files changed, 195 insertions(+), 19 deletions(-) diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp index 7d8beafb03e0..cc528801fe15 100644 --- a/SU2_CFD/include/limiters/computeLimiters_impl.hpp +++ b/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(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 < nPoint; ++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,21 @@ void computeLimiters_impl(CSolver* solver, for (size_t iVar = varBegin; iVar < varEnd; ++iVar) projMax[iVar] = projMin[iVar] = 0.0; + if (periodicProj != nullptr) + { + /*--- Start from the min/max over the edges of the periodic matches. ---*/ + + for (auto iVar = varBegin; iVar < varEnd; ++iVar) + { + const auto& periodicMin = (*periodicProj)(iPoint, iVar); + const auto& periodicMax = (*periodicProj)(iPoint, 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)) { @@ -224,7 +257,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) diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 7a7b6c726021..88e09d166dd1 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -146,6 +146,7 @@ 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 and max reconstruction increments over the edges of the rotational periodic matches of each point (for limiters). */ bool dynamic_grid; /*!< \brief Flag that determines whether the grid is dynamic (moving or deforming + grid velocities). */ @@ -4231,6 +4232,14 @@ 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] config - Definition of the particular problem. + * \return The matrix (nPoint x 2*nPrimVarGrad, min then max), nullptr if no periodic marker rotates the solution. + */ + su2activematrix* GetPeriodicProjections(const CConfig& config); + /*! * \brief Retrieve the solver name for output purposes. * \returns Name of the solver. diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 2f62f3241f36..e65f3f0ba0a4 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -262,7 +262,8 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, JCOUNT = nDim; break; case PERIODIC_LIM_PRIM_1: - COUNT_PER_POINT = nPrimVarGrad*2; + /*--- Min and max of the solution, min and max of the reconstruction increments. ---*/ + COUNT_PER_POINT = nPrimVarGrad*4; MPI_TYPE = COMM_TYPE::DOUBLE; ICOUNT = nPrimVarGrad; break; @@ -336,6 +337,31 @@ namespace PeriodicCommHelpers { break; } } + + /*--- Whether a periodic marker with these angles rotates the vector components of the solution. ---*/ + bool isRotation(const su2double* angles) { + return angles[0] != 0.0 || angles[1] != 0.0 || angles[2] != 0.0; + } +} + +su2activematrix* CSolver::GetPeriodicProjections(const CConfig& config) { + + if (!rotate_periodic) return nullptr; + + bool rotation = false; + for (auto iMarker = 0u; iMarker < config.GetnMarker_All(); iMarker++) { + if (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) continue; + rotation |= PeriodicCommHelpers::isRotation(config.GetPeriodicRotAngles(config.GetMarker_All_TagBound(iMarker))); + } + if (!rotation) return nullptr; + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { + if (PeriodicProj.rows() != nPoint) PeriodicProj.resize(nPoint, 2*nPrimVarGrad) = su2double(0.0); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + return &PeriodicProj; } void CSolver::InitiatePeriodicComms(CGeometry *geometry, @@ -391,6 +417,25 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, else GeometryToolbox::Rotate(rotMatrix3D, origin, direction, rotated); }; + /*--- Rotates, in place, the range [vMin, vMax] of each component of a vector: each term + of the product with the rotation matrix takes its smallest and its largest value. ---*/ + auto RotateBox = [&](su2double* vMin, su2double* vMax) { + su2double rotMin[3] = {0.0}, rotMax[3] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) { + for (auto jDim = 0u; jDim < nDim; jDim++) { + const su2double rotMatrix_ij = (nDim == 2) ? rotMatrix2D[iDim][jDim] : rotMatrix3D[iDim][jDim]; + const su2double fromMin = rotMatrix_ij * vMin[jDim]; + const su2double fromMax = rotMatrix_ij * vMax[jDim]; + rotMin[iDim] += min(fromMin, fromMax); + rotMax[iDim] += max(fromMin, fromMax); + } + } + for (auto iDim = 0u; iDim < nDim; iDim++) { + vMin[iDim] = rotMin[iDim]; + vMax[iDim] = rotMax[iDim]; + } + }; + string Marker_Tag; /*--- Set the size of the data packet and type depending on quantity. ---*/ @@ -421,6 +466,36 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, auto& limiter = PeriodicCommHelpers::selectLimiter(base_nodes, commType); auto& field = PeriodicCommHelpers::selectField(base_nodes, commType); + /*--- Updates the range [Sol_Min, Sol_Max] of each variable with the values of a neighbor, + the velocity of the neighbor (variables 1 to nDim) is rotated first if requested. ---*/ + auto UpdateMinMax = [&](const su2double* values, bool rotateVelocity) { + if (rotateVelocity) { + for (auto iVar = 0u; iVar < ICOUNT; iVar++) rotPrim_j[iVar] = values[iVar]; + Rotate(zeros, &values[1], &rotPrim_j[1]); + values = rotPrim_j; + } + for (auto iVar = 0u; iVar < ICOUNT; iVar++) { + Sol_Min[iVar] = min(Sol_Min[iVar], values[iVar]); + Sol_Max[iVar] = max(Sol_Max[iVar], values[iVar]); + } + }; + + /*--- Reconstruction increments of all variables from a point to the middle + of the edge to one of its neighbors, computed as in computeLimiters_impl. ---*/ + auto ReconstructionIncrements = [&](unsigned long point_i, unsigned long point_j, su2double* increments) { + const su2double kappa = config->GetMUSCL_Kappa_Flow(); + const auto* coord_i = geometry->nodes->GetCoord(point_i); + const auto* coord_j = geometry->nodes->GetCoord(point_j); + + for (auto iVar = 0u; iVar < ICOUNT; iVar++) { + su2double proj = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) + proj += 0.5 * (coord_j[iDim] - coord_i[iDim]) * gradient(point_i, iVar, iDim); + const su2double cent = 0.5 * (field(point_j, iVar) - field(point_i, iVar)); + increments[iVar] = LimiterHelpers<>::umusclProjection(proj, cent, kappa); + } + }; + /*--- Load the specified quantity from the solver into the generic communication buffer in the geometry class. ---*/ @@ -479,6 +554,10 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, GeometryToolbox::RotationMatrix(Theta, Phi, Psi, rotMatrix3D); } + /*--- Whether the vector components of the solution are rotated for this marker. ---*/ + + const bool rotation = rotate_periodic && PeriodicCommHelpers::isRotation(angles); + /*--- Compute the offset in the recv buffer for this point. ---*/ buf_offset = (msg_offset + iSend)*COUNT_PER_POINT; @@ -945,11 +1024,15 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, Sol_Max[iVar] = base_nodes->GetSolution_Max()(iPoint, iVar); } + /*--- The velocity of each neighbour is rotated before taking the min and max of its + components (the rotated min and max vectors do not bound the rotated components). + The min/max stored so far (the value of the point, or what was accumulated from + other periodic pairs) are a range for each component, which is rotated as a box. ---*/ + + if (rotation) RotateBox(&Sol_Min[1], &Sol_Max[1]); + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { - for (iVar = 0; iVar < ICOUNT; iVar++) { - Sol_Min[iVar] = min(Sol_Min[iVar], field(jPoint, iVar)); - Sol_Max[iVar] = max(Sol_Max[iVar], field(jPoint, iVar)); - } + UpdateMinMax(field[jPoint], rotation); } for (iVar = 0; iVar < ICOUNT; iVar++) { @@ -957,13 +1040,36 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, bufDSend[buf_offset+ICOUNT+iVar] = Sol_Max[iVar]; } - /*--- Rotate the momentum components of the min/max. ---*/ - - if (rotate_periodic) { + /*--- Preserve the original arithmetic for translational periodicity, including its AD recording. ---*/ + if (rotate_periodic && !rotation) { Rotate(zeros, &Sol_Min[1], &bufDSend[buf_offset+1]); Rotate(zeros, &Sol_Max[1], &bufDSend[buf_offset+ICOUNT+1]); } + /*--- The limiters of the velocity components cannot be compared across a rotation + (second phase). Instead, we also send the min and max, over "our" edges, of the + reconstruction increments (computed as in computeLimiters_impl) with the velocity + part rotated, the periodic match then computes its limiters with the complete + stencil. Without rotation these are zero, i.e. they have no effect. ---*/ + + if (commType == PERIODIC_LIM_PRIM_1) { + + for (auto iVar = 0u; iVar < ICOUNT; iVar++) + Sol_Min[iVar] = Sol_Max[iVar] = 0.0; + + if (rotation) { + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + ReconstructionIncrements(iPoint, jPoint, rotPrim_i); + UpdateMinMax(rotPrim_i, true); + } + } + + for (auto iVar = 0u; iVar < ICOUNT; iVar++) { + bufDSend[buf_offset+2*ICOUNT+iVar] = Sol_Min[iVar]; + bufDSend[buf_offset+3*ICOUNT+iVar] = Sol_Max[iVar]; + } + } + break; case PERIODIC_LIM_PRIM_2: @@ -977,7 +1083,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, bufDSend[buf_offset+iVar] = limiter(iPoint, iVar); } - if (rotate_periodic) { + /*--- The limiters of the velocity components cannot be compared across a rotation, + they were computed with the complete stencil (see the first phase), we send + a value that is never the minimum. ---*/ + + if (rotation) { + for (auto iDim = 0u; iDim < nDim; iDim++) + bufDSend[buf_offset+1+iDim] = std::numeric_limits::max(); + } else if (rotate_periodic) { Rotate(zeros, &limiter(iPoint,1), &bufDSend[buf_offset+1]); } @@ -1298,6 +1411,18 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, base_nodes->GetSolution_Max()(iPoint, iVar) = Solution_Max; } + /*--- Min/max reconstruction increments over the edges of the periodic match, + used to start the search over "our" edges (only with rotation). ---*/ + + if ((commType == PERIODIC_LIM_PRIM_1) && !PeriodicProj.empty()) { + for (auto iVar = 0u; iVar < ICOUNT; iVar++) { + PeriodicProj(iPoint, iVar) = min(PeriodicProj(iPoint, iVar), + bufDRecv[buf_offset+2*ICOUNT+iVar]); + PeriodicProj(iPoint, ICOUNT+iVar) = max(PeriodicProj(iPoint, ICOUNT+iVar), + bufDRecv[buf_offset+3*ICOUNT+iVar]); + } + } + break; case PERIODIC_LIM_PRIM_2: diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index f63748ba9283..fe259299d955 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -158,7 +158,7 @@ def main(): periodic2d.cfg_dir = "navierstokes/periodic2D" periodic2d.cfg_file = "config.cfg" periodic2d.test_iter = 1400 - periodic2d.test_vals = [-10.817607, -8.363541, -8.287457, -5.334100, -1.088412, -2945.200000] + periodic2d.test_vals = [-10.314545, -8.056213, -7.708669, -4.833054, -1.059806, -2953.400000] test_list.append(periodic2d) ########################## @@ -602,7 +602,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.907561, -12.214137, -19.151426, -13.451780, -19.085218, -13.450796, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.925871, -12.203300, -19.179088, -13.470092, -19.030151, -13.393012, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -627,7 +627,7 @@ def main(): multi_interface.cfg_dir = "turbomachinery/multi_interface" multi_interface.cfg_file = "multi_interface_rst.cfg" multi_interface.test_iter = 5 - multi_interface.test_vals = [-8.632240, -8.894740, -9.348706] + multi_interface.test_vals = [-8.634571, -8.895558, -9.348754] multi_interface.test_vals_aarch64 = [-8.632229, -8.894737, -9.348730] test_list.append(multi_interface) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 5bf6f7fb81f3..923a90bda73b 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -422,21 +422,29 @@ def main(): test_list.append(poiseuille_profile) # 2D rotational periodic sector without limiter (Jacobian of the periodic points) - periodic2d_no_limiter = TestCase('periodic2d_no_limiter') - periodic2d_no_limiter.cfg_dir = "navierstokes/periodic2D" - periodic2d_no_limiter.cfg_file = "no_limiter.cfg" + periodic2d_no_limiter = TestCase('periodic2d_no_limiter') + periodic2d_no_limiter.cfg_dir = "navierstokes/periodic2D" + periodic2d_no_limiter.cfg_file = "no_limiter.cfg" periodic2d_no_limiter.test_iter = 100 periodic2d_no_limiter.test_vals = [-3.216881, -0.582813, -0.629049, 2.266141, -1.056877, -811.470000] test_list.append(periodic2d_no_limiter) # 2D rotational periodic sector with multigrid - periodic2d_multigrid = TestCase('periodic2d_multigrid') - periodic2d_multigrid.cfg_dir = "navierstokes/periodic2D" - periodic2d_multigrid.cfg_file = "multigrid.cfg" + periodic2d_multigrid = TestCase('periodic2d_multigrid') + periodic2d_multigrid.cfg_dir = "navierstokes/periodic2D" + periodic2d_multigrid.cfg_file = "multigrid.cfg" periodic2d_multigrid.test_iter = 300 periodic2d_multigrid.test_vals = [-4.422778, -1.615251, -1.453736, 1.047350, -1.410225, -2060.600000] test_list.append(periodic2d_multigrid) + # 2D rotational periodic sector with limiter + periodic2d = TestCase('periodic2d') + periodic2d.cfg_dir = "navierstokes/periodic2D" + periodic2d.cfg_file = "config.cfg" + periodic2d.test_iter = 100 + periodic2d.test_vals = [-3.266248, -0.617784, -0.620252, 2.216476, -1.007761, -1035.000000] + test_list.append(periodic2d) + ########################## ### Compressible RANS ### ########################## From 3f9ac2ab49e1cefb7bdd0e404f1508ba475015c5 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 15:00:58 +0200 Subject: [PATCH 04/11] Refresh Jones restart references from periodic CI --- TestCases/hybrid_regression.py | 2 +- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index fe259299d955..cc24dd1f2601 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -602,7 +602,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.925871, -12.203300, -19.179088, -13.470092, -19.030151, -13.393012, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.959252, -12.213981, -19.314222, -13.600248, -19.076384, -13.444335, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 923a90bda73b..87b0f125ff51 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1276,7 +1276,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.910586, -12.203941, -19.198816, -13.489664, -19.033845, -13.400589, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.905382, -12.212100, -19.212143, -13.508268, -19.078783, -13.446389, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index e1a8c131e34e..e7c0261da35f 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -958,7 +958,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.916186, -12.202434, -19.225637, -13.511423, -19.026840, -13.394199, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.905642, -12.212181, -19.212554, -13.508053, -19.069991, -13.439494, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage From 9de9fd0b8298cc17a6e0fe4df8c2a4f74f8a3159 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 16:58:35 +0200 Subject: [PATCH 05/11] Avoid shadowing periodic communication loop indices --- SU2_CFD/src/solvers/CSolver.cpp | 56 ++++++++++++++++----------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index e65f3f0ba0a4..a8d90309eb3b 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -421,18 +421,18 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, of the product with the rotation matrix takes its smallest and its largest value. ---*/ auto RotateBox = [&](su2double* vMin, su2double* vMax) { su2double rotMin[3] = {0.0}, rotMax[3] = {0.0}; - for (auto iDim = 0u; iDim < nDim; iDim++) { + for (auto iCoordinate = 0u; iCoordinate < nDim; iCoordinate++) { for (auto jDim = 0u; jDim < nDim; jDim++) { - const su2double rotMatrix_ij = (nDim == 2) ? rotMatrix2D[iDim][jDim] : rotMatrix3D[iDim][jDim]; + const su2double rotMatrix_ij = (nDim == 2) ? rotMatrix2D[iCoordinate][jDim] : rotMatrix3D[iCoordinate][jDim]; const su2double fromMin = rotMatrix_ij * vMin[jDim]; const su2double fromMax = rotMatrix_ij * vMax[jDim]; - rotMin[iDim] += min(fromMin, fromMax); - rotMax[iDim] += max(fromMin, fromMax); + rotMin[iCoordinate] += min(fromMin, fromMax); + rotMax[iCoordinate] += max(fromMin, fromMax); } } - for (auto iDim = 0u; iDim < nDim; iDim++) { - vMin[iDim] = rotMin[iDim]; - vMax[iDim] = rotMax[iDim]; + for (auto iCoordinate = 0u; iCoordinate < nDim; iCoordinate++) { + vMin[iCoordinate] = rotMin[iCoordinate]; + vMax[iCoordinate] = rotMax[iCoordinate]; } }; @@ -470,13 +470,13 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, the velocity of the neighbor (variables 1 to nDim) is rotated first if requested. ---*/ auto UpdateMinMax = [&](const su2double* values, bool rotateVelocity) { if (rotateVelocity) { - for (auto iVar = 0u; iVar < ICOUNT; iVar++) rotPrim_j[iVar] = values[iVar]; + for (auto iField = 0u; iField < ICOUNT; iField++) rotPrim_j[iField] = values[iField]; Rotate(zeros, &values[1], &rotPrim_j[1]); values = rotPrim_j; } - for (auto iVar = 0u; iVar < ICOUNT; iVar++) { - Sol_Min[iVar] = min(Sol_Min[iVar], values[iVar]); - Sol_Max[iVar] = max(Sol_Max[iVar], values[iVar]); + for (auto iField = 0u; iField < ICOUNT; iField++) { + Sol_Min[iField] = min(Sol_Min[iField], values[iField]); + Sol_Max[iField] = max(Sol_Max[iField], values[iField]); } }; @@ -487,12 +487,12 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, const auto* coord_i = geometry->nodes->GetCoord(point_i); const auto* coord_j = geometry->nodes->GetCoord(point_j); - for (auto iVar = 0u; iVar < ICOUNT; iVar++) { + for (auto iField = 0u; iField < ICOUNT; iField++) { su2double proj = 0.0; - for (auto iDim = 0u; iDim < nDim; iDim++) - proj += 0.5 * (coord_j[iDim] - coord_i[iDim]) * gradient(point_i, iVar, iDim); - const su2double cent = 0.5 * (field(point_j, iVar) - field(point_i, iVar)); - increments[iVar] = LimiterHelpers<>::umusclProjection(proj, cent, kappa); + for (auto iCoordinate = 0u; iCoordinate < nDim; iCoordinate++) + proj += 0.5 * (coord_j[iCoordinate] - coord_i[iCoordinate]) * gradient(point_i, iField, iCoordinate); + const su2double cent = 0.5 * (field(point_j, iField) - field(point_i, iField)); + increments[iField] = LimiterHelpers<>::umusclProjection(proj, cent, kappa); } }; @@ -1054,8 +1054,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, if (commType == PERIODIC_LIM_PRIM_1) { - for (auto iVar = 0u; iVar < ICOUNT; iVar++) - Sol_Min[iVar] = Sol_Max[iVar] = 0.0; + for (auto iField = 0u; iField < ICOUNT; iField++) + Sol_Min[iField] = Sol_Max[iField] = 0.0; if (rotation) { for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { @@ -1064,9 +1064,9 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, } } - for (auto iVar = 0u; iVar < ICOUNT; iVar++) { - bufDSend[buf_offset+2*ICOUNT+iVar] = Sol_Min[iVar]; - bufDSend[buf_offset+3*ICOUNT+iVar] = Sol_Max[iVar]; + for (auto iField = 0u; iField < ICOUNT; iField++) { + bufDSend[buf_offset+2*ICOUNT+iField] = Sol_Min[iField]; + bufDSend[buf_offset+3*ICOUNT+iField] = Sol_Max[iField]; } } @@ -1088,8 +1088,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, a value that is never the minimum. ---*/ if (rotation) { - for (auto iDim = 0u; iDim < nDim; iDim++) - bufDSend[buf_offset+1+iDim] = std::numeric_limits::max(); + for (auto iCoordinate = 0u; iCoordinate < nDim; iCoordinate++) + bufDSend[buf_offset+1+iCoordinate] = std::numeric_limits::max(); } else if (rotate_periodic) { Rotate(zeros, &limiter(iPoint,1), &bufDSend[buf_offset+1]); } @@ -1415,11 +1415,11 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, used to start the search over "our" edges (only with rotation). ---*/ if ((commType == PERIODIC_LIM_PRIM_1) && !PeriodicProj.empty()) { - for (auto iVar = 0u; iVar < ICOUNT; iVar++) { - PeriodicProj(iPoint, iVar) = min(PeriodicProj(iPoint, iVar), - bufDRecv[buf_offset+2*ICOUNT+iVar]); - PeriodicProj(iPoint, ICOUNT+iVar) = max(PeriodicProj(iPoint, ICOUNT+iVar), - bufDRecv[buf_offset+3*ICOUNT+iVar]); + for (auto iField = 0u; iField < ICOUNT; iField++) { + PeriodicProj(iPoint, iField) = min(PeriodicProj(iPoint, iField), + bufDRecv[buf_offset+2*ICOUNT+iField]); + PeriodicProj(iPoint, ICOUNT+iField) = max(PeriodicProj(iPoint, ICOUNT+iField), + bufDRecv[buf_offset+3*ICOUNT+iField]); } } From c5ebebec27da40aec22afd64814a753430430281 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 8 Oct 2026 20:28:52 +0200 Subject: [PATCH 06/11] Address rotational periodic limiter review feedback --- Common/include/toolboxes/geometry_toolbox.hpp | 25 ++++ SU2_CFD/include/limiters/CLimiterDetails.hpp | 12 ++ .../include/limiters/computeLimiters_impl.hpp | 40 ++---- SU2_CFD/include/solvers/CSolver.hpp | 16 ++- SU2_CFD/src/solvers/CSolver.cpp | 49 +++---- .../navierstokes/periodic2D/multigrid.cfg | 10 +- .../navierstokes/periodic2D/no_limiter.cfg | 10 +- UnitTests/SU2_CFD/periodic_limiters.cpp | 129 ++++++++++++++++++ UnitTests/meson.build | 1 + 9 files changed, 230 insertions(+), 62 deletions(-) create mode 100644 UnitTests/SU2_CFD/periodic_limiters.cpp diff --git a/Common/include/toolboxes/geometry_toolbox.hpp b/Common/include/toolboxes/geometry_toolbox.hpp index 1087f940a36f..05995fbc84c2 100644 --- a/Common/include/toolboxes/geometry_toolbox.hpp +++ b/Common/include/toolboxes/geometry_toolbox.hpp @@ -27,6 +27,7 @@ #pragma once #include +#include namespace GeometryToolbox { /// \addtogroup GeometryToolbox @@ -215,6 +216,30 @@ 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 +inline bool HasRotation(const Scalar* angles) { + return angles[0] != 0.0 || angles[1] != 0.0 || angles[2] != 0.0; +} + +/*! \brief Rotate component bounds in place, enclosing the rotated box. */ +template +inline void RotateBox(const Scalar R[][nDim], Scalar* vMin, Scalar* vMax) { + 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] += std::min(fromMin, fromMax); + rotMax[iDim] += std::max(fromMin, fromMax); + } + } + for (int iDim = 0; iDim < nDim; ++iDim) { + vMin[iDim] = rotMin[iDim]; + vMax[iDim] = rotMax[iDim]; + } +} + /*! \brief Tangent projection */ template inline void TangentProjection(Int nDim, const Mat& tensor, const Scalar* vector, Scalar* proj) { diff --git a/SU2_CFD/include/limiters/CLimiterDetails.hpp b/SU2_CFD/include/limiters/CLimiterDetails.hpp index 2fe298261006..ca2577c78c19 100644 --- a/SU2_CFD/include/limiters/CLimiterDetails.hpp +++ b/SU2_CFD/include/limiters/CLimiterDetails.hpp @@ -72,6 +72,18 @@ struct LimiterHelpers { FORCEINLINE static Type epsilon() {return std::numeric_limits::epsilon();} + /*! \brief MUSCL reconstruction increment from a point to the middle of an edge. */ + template + FORCEINLINE static Type reconstructionIncrement(Int nDim, const Type* coord_i, const Type* coord_j, + 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 += 0.5 * (coord_j[iDim] - coord_i[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) { /*-------------------------------------------------------------------*/ diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp index cc528801fe15..2559649168cf 100644 --- a/SU2_CFD/include/limiters/computeLimiters_impl.hpp +++ b/SU2_CFD/include/limiters/computeLimiters_impl.hpp @@ -116,7 +116,7 @@ void computeLimiters_impl(CSolver* solver, su2activematrix* periodicProj = nullptr; if (periodic && (kindPeriodicComm1 == PERIODIC_LIM_PRIM_1)) - periodicProj = solver->GetPeriodicProjections(config); + periodicProj = solver->GetPeriodicProjections(geometry, config); /*--- Initialize all min/max field values if we have * periodic comms. otherwise do it inside main loop. ---*/ @@ -126,7 +126,7 @@ void computeLimiters_impl(CSolver* solver, if (periodicProj != nullptr) { SU2_OMP_FOR_STAT(chunkSize) - for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) + 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 @@ -183,18 +183,20 @@ void computeLimiters_impl(CSolver* solver, for (size_t iVar = varBegin; iVar < varEnd; ++iVar) projMax[iVar] = projMin[iVar] = 0.0; - if (periodicProj != nullptr) + if (periodicProj != nullptr && nodes->GetPeriodicBoundary(iPoint)) { /*--- Start from the min/max over the edges of the periodic matches. ---*/ - for (auto iVar = varBegin; iVar < varEnd; ++iVar) - { - const auto& periodicMin = (*periodicProj)(iPoint, iVar); - const auto& periodicMax = (*periodicProj)(iPoint, periodicProj->cols()/2 + iVar); - AD::SetPreaccIn(periodicMin); - AD::SetPreaccIn(periodicMax); - projMin[iVar] = periodicMin; - projMax[iVar] = periodicMax; + 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; + } } } @@ -205,25 +207,13 @@ void computeLimiters_impl(CSolver* solver, const auto coord_j = geometry.nodes->GetCoord(jPoint); AD::SetPreaccIn(coord_j, nDim); - /*--- Distance vector from iPoint to face (middle of the edge). ---*/ - - su2double dist_ij[nDim] = {0.0}; - - 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, coord_i, coord_j, + gradient[iPoint][iVar], field(iPoint,iVar), field(jPoint,iVar), umusclKappa); projMax[iVar] = max(projMax[iVar], proj); projMin[iVar] = min(projMin[iVar], proj); diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 88e09d166dd1..76d883f2690b 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -37,6 +37,7 @@ #include #include #include +#include #include #include @@ -146,7 +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 and max reconstruction increments over the edges of the rotational periodic matches of each point (for limiters). */ + su2activematrix PeriodicProj; /*!< \brief Min/max reconstruction increments at periodic receive points (for limiters). */ + std::unordered_map 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). */ @@ -4235,10 +4237,18 @@ class CSolver { /*! * \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 (nPoint x 2*nPrimVarGrad, min then max), nullptr if no periodic marker rotates the solution. + * \return The matrix (unique periodic receive points x 2*nPrimVarGrad, min then max), nullptr without rotation. */ - su2activematrix* GetPeriodicProjections(const CConfig& config); + 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. diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index a8d90309eb3b..407a287612da 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -337,27 +337,26 @@ namespace PeriodicCommHelpers { break; } } - - /*--- Whether a periodic marker with these angles rotates the vector components of the solution. ---*/ - bool isRotation(const su2double* angles) { - return angles[0] != 0.0 || angles[1] != 0.0 || angles[2] != 0.0; - } } -su2activematrix* CSolver::GetPeriodicProjections(const CConfig& config) { +su2activematrix* CSolver::GetPeriodicProjections(const CGeometry& geometry, const CConfig& config) { if (!rotate_periodic) return nullptr; bool rotation = false; for (auto iMarker = 0u; iMarker < config.GetnMarker_All(); iMarker++) { if (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) continue; - rotation |= PeriodicCommHelpers::isRotation(config.GetPeriodicRotAngles(config.GetMarker_All_TagBound(iMarker))); + rotation |= GeometryToolbox::HasRotation(config.GetPeriodicRotAngles(config.GetMarker_All_TagBound(iMarker))); } if (!rotation) return nullptr; BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - if (PeriodicProj.rows() != nPoint) PeriodicProj.resize(nPoint, 2*nPrimVarGrad) = su2double(0.0); + if (PeriodicProj.empty() && geometry.nPeriodicRecv > 0) { + for (auto iRecv = 0; iRecv < geometry.nPoint_PeriodicRecv[geometry.nPeriodicRecv]; ++iRecv) + PeriodicProjIndex.emplace(geometry.Local_Point_PeriodicRecv[iRecv], PeriodicProjIndex.size()); + PeriodicProj.resize(PeriodicProjIndex.size(), 2*nPrimVarGrad); + } } END_SU2_OMP_SAFE_GLOBAL_ACCESS @@ -420,20 +419,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Rotates, in place, the range [vMin, vMax] of each component of a vector: each term of the product with the rotation matrix takes its smallest and its largest value. ---*/ auto RotateBox = [&](su2double* vMin, su2double* vMax) { - su2double rotMin[3] = {0.0}, rotMax[3] = {0.0}; - for (auto iCoordinate = 0u; iCoordinate < nDim; iCoordinate++) { - for (auto jDim = 0u; jDim < nDim; jDim++) { - const su2double rotMatrix_ij = (nDim == 2) ? rotMatrix2D[iCoordinate][jDim] : rotMatrix3D[iCoordinate][jDim]; - const su2double fromMin = rotMatrix_ij * vMin[jDim]; - const su2double fromMax = rotMatrix_ij * vMax[jDim]; - rotMin[iCoordinate] += min(fromMin, fromMax); - rotMax[iCoordinate] += max(fromMin, fromMax); - } - } - for (auto iCoordinate = 0u; iCoordinate < nDim; iCoordinate++) { - vMin[iCoordinate] = rotMin[iCoordinate]; - vMax[iCoordinate] = rotMax[iCoordinate]; - } + if (nDim == 2) GeometryToolbox::RotateBox(rotMatrix2D, vMin, vMax); + else GeometryToolbox::RotateBox(rotMatrix3D, vMin, vMax); }; string Marker_Tag; @@ -488,11 +475,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, const auto* coord_j = geometry->nodes->GetCoord(point_j); for (auto iField = 0u; iField < ICOUNT; iField++) { - su2double proj = 0.0; - for (auto iCoordinate = 0u; iCoordinate < nDim; iCoordinate++) - proj += 0.5 * (coord_j[iCoordinate] - coord_i[iCoordinate]) * gradient(point_i, iField, iCoordinate); - const su2double cent = 0.5 * (field(point_j, iField) - field(point_i, iField)); - increments[iField] = LimiterHelpers<>::umusclProjection(proj, cent, kappa); + increments[iField] = LimiterHelpers<>::reconstructionIncrement(nDim, coord_i, coord_j, + gradient[point_i][iField], field(point_i, iField), field(point_j, iField), kappa); } }; @@ -556,7 +540,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Whether the vector components of the solution are rotated for this marker. ---*/ - const bool rotation = rotate_periodic && PeriodicCommHelpers::isRotation(angles); + const bool rotation = rotate_periodic && GeometryToolbox::HasRotation(angles); /*--- Compute the offset in the recv buffer for this point. ---*/ @@ -1415,11 +1399,12 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, used to start the search over "our" edges (only with rotation). ---*/ if ((commType == PERIODIC_LIM_PRIM_1) && !PeriodicProj.empty()) { + auto* projections = GetPeriodicProjection(iPoint); + assert(projections != nullptr); for (auto iField = 0u; iField < ICOUNT; iField++) { - PeriodicProj(iPoint, iField) = min(PeriodicProj(iPoint, iField), - bufDRecv[buf_offset+2*ICOUNT+iField]); - PeriodicProj(iPoint, ICOUNT+iField) = max(PeriodicProj(iPoint, ICOUNT+iField), - bufDRecv[buf_offset+3*ICOUNT+iField]); + projections[iField] = min(projections[iField], bufDRecv[buf_offset+2*ICOUNT+iField]); + projections[ICOUNT+iField] = max(projections[ICOUNT+iField], + bufDRecv[buf_offset+3*ICOUNT+iField]); } } diff --git a/TestCases/navierstokes/periodic2D/multigrid.cfg b/TestCases/navierstokes/periodic2D/multigrid.cfg index 72c8e9933439..d7a7c67bbb31 100644 --- a/TestCases/navierstokes/periodic2D/multigrid.cfg +++ b/TestCases/navierstokes/periodic2D/multigrid.cfg @@ -1,5 +1,13 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Same as config.cfg with multigrid, without slope limiter, lower CFL. +% % +% SU2 configuration file % +% Case description: Rotational periodic sector with multigrid, % +% without slope limiter, at lower CFL % +% Author: SU2 Contributors % +% Date: Oct 2026 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SOLVER= NAVIER_STOKES KIND_TURB_MODEL= NONE diff --git a/TestCases/navierstokes/periodic2D/no_limiter.cfg b/TestCases/navierstokes/periodic2D/no_limiter.cfg index 1c0a9320b94b..2ac3b09dfddb 100644 --- a/TestCases/navierstokes/periodic2D/no_limiter.cfg +++ b/TestCases/navierstokes/periodic2D/no_limiter.cfg @@ -1,5 +1,13 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Same as config.cfg without slope limiter (tests the Jacobian of the periodic points). +% % +% SU2 configuration file % +% Case description: Rotational periodic sector without slope limiter, % +% testing the periodic Jacobian % +% Author: SU2 Contributors % +% Date: Oct 2026 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SOLVER= NAVIER_STOKES KIND_TURB_MODEL= NONE diff --git a/UnitTests/SU2_CFD/periodic_limiters.cpp b/UnitTests/SU2_CFD/periodic_limiters.cpp new file mode 100644 index 000000000000..9dbd686c7d08 --- /dev/null +++ b/UnitTests/SU2_CFD/periodic_limiters.cpp @@ -0,0 +1,129 @@ +/*! + * \file periodic_limiters.cpp + * \brief Tests for periodic limiter storage and reconstruction helpers. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include "../../SU2_CFD/include/solvers/CSolver.hpp" + +TEST_CASE("Rotation of component bounds encloses every corner", "[PeriodicLimiter]") { + const su2double noAngles[3] = {0.0, 0.0, 0.0}, negativeAngle[3] = {0.0, -0.5, 0.0}; + CHECK_FALSE(GeometryToolbox::HasRotation(noAngles)); + CHECK(GeometryToolbox::HasRotation(negativeAngle)); + auto checkCorners = [](auto& rotation, auto& lower, auto& upper) { + constexpr int nDim = sizeof(lower) / sizeof(lower[0]); + su2double expectedMin[nDim], expectedMax[nDim]; + for (int iDim = 0; iDim < nDim; ++iDim) { + expectedMin[iDim] = std::numeric_limits::max(); + expectedMax[iDim] = -expectedMin[iDim]; + } + const su2double origin[nDim] = {0.0}; + for (int corner = 0; corner < (1 << nDim); ++corner) { + su2double point[nDim], rotated[nDim]; + for (int iDim = 0; iDim < nDim; ++iDim) point[iDim] = (corner & (1 << iDim)) ? upper[iDim] : lower[iDim]; + GeometryToolbox::Rotate(rotation, origin, point, rotated); + for (int iDim = 0; iDim < nDim; ++iDim) { + expectedMin[iDim] = std::min(expectedMin[iDim], rotated[iDim]); + expectedMax[iDim] = std::max(expectedMax[iDim], rotated[iDim]); + } + } + GeometryToolbox::RotateBox(rotation, lower, upper); + for (int iDim = 0; iDim < nDim; ++iDim) { + CHECK(lower[iDim] == Approx(expectedMin[iDim])); + CHECK(upper[iDim] == Approx(expectedMax[iDim])); + } + }; + + SECTION("2D oblique rotation") { + su2double rotation[2][2], lower[2] = {-2.0, 1.0}, upper[2] = {3.0, 4.0}; + GeometryToolbox::RotationMatrix(su2double(PI_NUMBER / 4.0), rotation); + checkCorners(rotation, lower, upper); + } + SECTION("3D rotation about all axes") { + su2double rotation[3][3], lower[3] = {-2.0, 1.0, -3.0}, upper[3] = {3.0, 4.0, 2.0}; + GeometryToolbox::RotationMatrix(su2double(0.3), su2double(-0.7), su2double(1.2), rotation); + checkCorners(rotation, lower, upper); + } +} + +TEST_CASE("Shared reconstruction increment has the MUSCL scaling", "[PeriodicLimiter]") { + const su2double coord_i[3] = {1.0, 2.0, 3.0}, coord_j[3] = {3.0, 0.0, 4.0}; + const su2double gradient[3] = {2.0, -1.0, 3.0}; + for (const auto kappa : {-1.0, 0.0, 0.5, 1.0}) { + const auto increment = LimiterHelpers<>::reconstructionIncrement(3, coord_i, coord_j, gradient, su2double(7.0), + su2double(15.0), su2double(kappa)); + CHECK(increment == Approx(4.5 - 0.5 * kappa)); + const auto linearIncrement = LimiterHelpers<>::reconstructionIncrement( + 3, coord_i, coord_j, gradient, su2double(7.0), su2double(16.0), su2double(kappa)); + CHECK(linearIncrement == Approx(4.5)); + } +} + +TEST_CASE("Periodic projection storage uses unique receive points", "[PeriodicLimiter]") { + std::stringstream options; + options << "SOLVER= EULER\n" + "MARKER_PERIODIC= (per1, per2, 0,0,0, 0,0,45, 0,0,0)\n"; + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + config.SetnMarker_All(1); + config.SetMarker_All_KindBC(0, PERIODIC_BOUNDARY); + config.SetMarker_All_TagBound(0, "per1"); + + struct ProjectionSolver : CSolver { + CVariable* GetBaseClassPointerToNodes() override { return nullptr; } + ProjectionSolver() { + nPoint = 1000000; + nVar = 4; + nPrimVarGrad = 4; + rotate_periodic = true; + } + } solver; + + CGeometry geometry; + geometry.nPeriodicRecv = 1; + geometry.nPoint_PeriodicRecv = new int[2]{0, 4}; + geometry.Local_Point_PeriodicRecv = new unsigned long[4]{12, 27, 12, 99999}; + + auto* storage = solver.GetPeriodicProjections(geometry, config); + REQUIRE(storage != nullptr); + CHECK(storage->rows() == 3); + CHECK(storage->cols() == 8); + CHECK(solver.GetPeriodicProjection(13) == nullptr); + REQUIRE(solver.GetPeriodicProjection(12) != nullptr); + REQUIRE(solver.GetPeriodicProjection(27) != nullptr); + REQUIRE(solver.GetPeriodicProjection(99999) != nullptr); + solver.GetPeriodicProjection(12)[0] = -2.0; + CHECK(solver.GetPeriodicProjection(27) != solver.GetPeriodicProjection(12)); + CHECK(solver.GetPeriodicProjections(geometry, config) == storage); + CHECK(solver.GetPeriodicProjection(12)[0] == -2.0); + + CGeometry noReceives; + ProjectionSolver solverWithoutReceives; + auto* emptyStorage = solverWithoutReceives.GetPeriodicProjections(noReceives, config); + REQUIRE(emptyStorage != nullptr); + CHECK(emptyStorage->empty()); + CHECK(solverWithoutReceives.GetPeriodicProjection(12) == nullptr); + + solver.SetRotatePeriodic(false); + CHECK(solver.GetPeriodicProjections(geometry, config) == nullptr); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index ba3c63afce92..1d1a97f3b827 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -19,6 +19,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'SU2_CFD/output/CNEMOCompOutput_tests.cpp', 'SU2_CFD/output/COutput_convergence_tests.cpp', 'SU2_CFD/gradients.cpp', + 'SU2_CFD/periodic_limiters.cpp', 'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp', 'SU2_CFD/nemo_viscous_assembly.cpp', 'SU2_CFD/windowing.cpp', From 609f0c2e4c4f7480ed2783a5f41b9a1bfe6c8760 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 8 Oct 2026 21:03:23 +0200 Subject: [PATCH 07/11] Exercise Aachen turbine multigrid in serial and MPI regressions --- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- .../turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg | 4 ++++ 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 87b0f125ff51..f2d4c669b5ad 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1268,7 +1268,7 @@ def main(): Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 Aachen_3D_restart.tol = 0.00001 - Aachen_3D_restart.test_vals = [-7.701423, -8.504851, -6.014951, -6.472219, -5.802012, -4.609768, -5.550659, -5.300718, -3.804222, -5.255983, -5.763064, -3.609604, -2.229253, -2.880518, -0.563484] + Aachen_3D_restart.test_vals = [-7.688484, -8.467124, -6.035067, -6.425070, -5.819109, -4.597817, -5.528414, -5.304377, -3.819864, -5.242170, -5.753869, -3.630202, -2.213242, -2.860178, -0.576864] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index e7c0261da35f..e3fd07978fa4 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -949,7 +949,7 @@ def main(): Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 - Aachen_3D_restart.test_vals = [-7.701424, -8.504850, -6.014951, -6.472217, -5.802013, -4.609762, -5.550659, -5.300718, -3.804222, -5.255982, -5.763064, -3.609600, -2.229254, -2.880518, -0.563485] + Aachen_3D_restart.test_vals = [-7.689116, -8.463749, -6.036107, -6.426604, -5.832700, -4.599027, -5.527178, -5.306021, -3.821792, -5.242926, -5.761004, -3.631603, -2.213929, -2.871564, -0.578545] Aachen_3D_restart.enabled_with_asan = False test_list.append(Aachen_3D_restart) diff --git a/TestCases/turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg b/TestCases/turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg index 0aebb8cd4d42..9e1c0f04cc45 100755 --- a/TestCases/turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg +++ b/TestCases/turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg @@ -228,6 +228,10 @@ LINEAR_SOLVER_ERROR= 1E-4 % Max number of iterations of the linear solver for the implicit formulation LINEAR_SOLVER_ITER= 15 % +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 1 +% % ----------------------- SLOPE LIMITER DEFINITION ----------------------------% % % Coefficient for the limiter From 639317d728ae5b0bd053d95b3528ab0dcd3e3803 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 8 Oct 2026 22:01:11 +0200 Subject: [PATCH 08/11] Preserve limiter displacement reuse and AD math dispatch --- Common/include/toolboxes/geometry_toolbox.hpp | 6 ++++-- SU2_CFD/include/limiters/CLimiterDetails.hpp | 10 ++++------ SU2_CFD/include/limiters/computeLimiters_impl.hpp | 8 +++++++- SU2_CFD/src/solvers/CSolver.cpp | 6 +++++- UnitTests/SU2_CFD/periodic_limiters.cpp | 12 ++++++++---- 5 files changed, 28 insertions(+), 14 deletions(-) diff --git a/Common/include/toolboxes/geometry_toolbox.hpp b/Common/include/toolboxes/geometry_toolbox.hpp index 05995fbc84c2..0e635197b321 100644 --- a/Common/include/toolboxes/geometry_toolbox.hpp +++ b/Common/include/toolboxes/geometry_toolbox.hpp @@ -225,13 +225,15 @@ inline bool HasRotation(const Scalar* angles) { /*! \brief Rotate component bounds in place, enclosing the rotated box. */ template 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] += std::min(fromMin, fromMax); - rotMax[iDim] += std::max(fromMin, fromMax); + rotMin[iDim] += min(fromMin, fromMax); + rotMax[iDim] += max(fromMin, fromMax); } } for (int iDim = 0; iDim < nDim; ++iDim) { diff --git a/SU2_CFD/include/limiters/CLimiterDetails.hpp b/SU2_CFD/include/limiters/CLimiterDetails.hpp index ca2577c78c19..554509f066fc 100644 --- a/SU2_CFD/include/limiters/CLimiterDetails.hpp +++ b/SU2_CFD/include/limiters/CLimiterDetails.hpp @@ -72,14 +72,12 @@ struct LimiterHelpers { FORCEINLINE static Type epsilon() {return std::numeric_limits::epsilon();} - /*! \brief MUSCL reconstruction increment from a point to the middle of an edge. */ + /*! \brief MUSCL reconstruction increment using the displacement to the middle of an edge. */ template - FORCEINLINE static Type reconstructionIncrement(Int nDim, const Type* coord_i, const Type* coord_j, - const Type* gradient, const Type& value_i, - const Type& value_j, const Type& kappa) { + 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 += 0.5 * (coord_j[iDim] - coord_i[iDim]) * gradient[iDim]; + 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); } diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp index 2559649168cf..3b19f241f3bc 100644 --- a/SU2_CFD/include/limiters/computeLimiters_impl.hpp +++ b/SU2_CFD/include/limiters/computeLimiters_impl.hpp @@ -207,12 +207,18 @@ void computeLimiters_impl(CSolver* solver, const auto coord_j = geometry.nodes->GetCoord(jPoint); AD::SetPreaccIn(coord_j, nDim); + /*--- Distance vector from iPoint to face (middle of the edge). ---*/ + + su2double dist_ij[nDim] = {0.0}; + 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) { AD::SetPreaccIn(field(jPoint,iVar)); - const su2double proj = LimiterHelpers<>::reconstructionIncrement(nDim, coord_i, coord_j, + const su2double proj = LimiterHelpers<>::reconstructionIncrement(nDim, dist_ij, gradient[iPoint][iVar], field(iPoint,iVar), field(jPoint,iVar), umusclKappa); projMax[iVar] = max(projMax[iVar], proj); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 407a287612da..fd02e9956918 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -474,8 +474,12 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, const auto* coord_i = geometry->nodes->GetCoord(point_i); const auto* coord_j = geometry->nodes->GetCoord(point_j); + su2double dist_ij[3] = {0.0}; + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + dist_ij[iCoordinate] = 0.5 * (coord_j[iCoordinate] - coord_i[iCoordinate]); + for (auto iField = 0u; iField < ICOUNT; iField++) { - increments[iField] = LimiterHelpers<>::reconstructionIncrement(nDim, coord_i, coord_j, + increments[iField] = LimiterHelpers<>::reconstructionIncrement(nDim, dist_ij, gradient[point_i][iField], field(point_i, iField), field(point_j, iField), kappa); } }; diff --git a/UnitTests/SU2_CFD/periodic_limiters.cpp b/UnitTests/SU2_CFD/periodic_limiters.cpp index 9dbd686c7d08..a7fc3222a0c9 100644 --- a/UnitTests/SU2_CFD/periodic_limiters.cpp +++ b/UnitTests/SU2_CFD/periodic_limiters.cpp @@ -68,14 +68,18 @@ TEST_CASE("Rotation of component bounds encloses every corner", "[PeriodicLimite } TEST_CASE("Shared reconstruction increment has the MUSCL scaling", "[PeriodicLimiter]") { - const su2double coord_i[3] = {1.0, 2.0, 3.0}, coord_j[3] = {3.0, 0.0, 4.0}; + const su2double halfEdge[3] = {1.0, -1.0, 0.5}; const su2double gradient[3] = {2.0, -1.0, 3.0}; for (const auto kappa : {-1.0, 0.0, 0.5, 1.0}) { - const auto increment = LimiterHelpers<>::reconstructionIncrement(3, coord_i, coord_j, gradient, su2double(7.0), + const auto increment = LimiterHelpers<>::reconstructionIncrement(3, halfEdge, gradient, su2double(7.0), su2double(15.0), su2double(kappa)); CHECK(increment == Approx(4.5 - 0.5 * kappa)); - const auto linearIncrement = LimiterHelpers<>::reconstructionIncrement( - 3, coord_i, coord_j, gradient, su2double(7.0), su2double(16.0), su2double(kappa)); + CHECK(LimiterHelpers<>::reconstructionIncrement(2, halfEdge, gradient, su2double(7.0), su2double(15.0), + su2double(kappa)) == Approx(3.0 + kappa)); + CHECK(LimiterHelpers<>::reconstructionIncrement(2, halfEdge, gradient, su2double(7.0), su2double(13.0), + su2double(kappa)) == Approx(3.0)); + const auto linearIncrement = LimiterHelpers<>::reconstructionIncrement(3, halfEdge, gradient, su2double(7.0), + su2double(16.0), su2double(kappa)); CHECK(linearIncrement == Approx(4.5)); } } From 492d2ebdc42888527b70f58f52ca30b6ab1bce62 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 8 Oct 2026 22:01:11 +0200 Subject: [PATCH 09/11] Use Aachen multigrid references from serial and MPI CI --- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index f2d4c669b5ad..89c6eb9a7c72 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1268,7 +1268,7 @@ def main(): Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 Aachen_3D_restart.tol = 0.00001 - Aachen_3D_restart.test_vals = [-7.688484, -8.467124, -6.035067, -6.425070, -5.819109, -4.597817, -5.528414, -5.304377, -3.819864, -5.242170, -5.753869, -3.630202, -2.213242, -2.860178, -0.576864] + Aachen_3D_restart.test_vals = [-7.688483, -8.466346, -6.035067, -6.425967, -5.820125, -4.597817, -5.528297, -5.301466, -3.819862, -5.242356, -5.738805, -3.630202, -2.213240, -2.857919, -0.576864] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index e3fd07978fa4..bd3cd676ca33 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -949,7 +949,7 @@ def main(): Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 - Aachen_3D_restart.test_vals = [-7.689116, -8.463749, -6.036107, -6.426604, -5.832700, -4.599027, -5.527178, -5.306021, -3.821792, -5.242926, -5.761004, -3.631603, -2.213929, -2.871564, -0.578545] + Aachen_3D_restart.test_vals = [-7.689117, -8.476204, -6.036107, -6.427402, -5.831576, -4.599024, -5.530512, -5.305687, -3.821789, -5.243255, -5.754831, -3.631603, -2.213934, -2.895728, -0.578545] Aachen_3D_restart.enabled_with_asan = False test_list.append(Aachen_3D_restart) From 6280ff1741ecb49cdb66c5ac82aa328b1988d44b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 8 Oct 2026 22:36:26 +0200 Subject: [PATCH 10/11] Refresh parallel AD stator reference after limiter refactor --- TestCases/parallel_regression_AD.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 26a32c3e2762..b80601e77efe 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -251,7 +251,7 @@ def main(): discadj_trans_stator.cfg_dir = "disc_adj_turbomachinery/transonic_stator_2D" discadj_trans_stator.cfg_file = "transonic_stator.cfg" discadj_trans_stator.test_iter = 79 - discadj_trans_stator.test_vals = [79.000000, -7.555663, -10.335501, -10.356934, -13.629559] + discadj_trans_stator.test_vals = [79.000000, -7.555647, -10.335486, -10.356919, -13.629543] discadj_trans_stator.test_vals_aarch64 = [79.000000, -7.555647, -10.335486, -10.356919, -13.629543] test_list.append(discadj_trans_stator) From 12d95ac76d59cacec1497020ea5dc5ff59bbaf45 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 9 Oct 2026 20:39:38 +0200 Subject: [PATCH 11/11] Document exact periodic angle checks and simplify bound copies --- Common/include/toolboxes/geometry_toolbox.hpp | 7 +++---- UnitTests/SU2_CFD/periodic_limiters.cpp | 19 +++++++++++++++++-- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/Common/include/toolboxes/geometry_toolbox.hpp b/Common/include/toolboxes/geometry_toolbox.hpp index 0e635197b321..e621ea65b7b2 100644 --- a/Common/include/toolboxes/geometry_toolbox.hpp +++ b/Common/include/toolboxes/geometry_toolbox.hpp @@ -219,6 +219,7 @@ 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 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; } @@ -236,10 +237,8 @@ inline void RotateBox(const Scalar R[][nDim], Scalar* vMin, Scalar* vMax) { rotMax[iDim] += max(fromMin, fromMax); } } - for (int iDim = 0; iDim < nDim; ++iDim) { - vMin[iDim] = rotMin[iDim]; - vMax[iDim] = rotMax[iDim]; - } + std::copy_n(rotMin, nDim, vMin); + std::copy_n(rotMax, nDim, vMax); } /*! \brief Tangent projection */ diff --git a/UnitTests/SU2_CFD/periodic_limiters.cpp b/UnitTests/SU2_CFD/periodic_limiters.cpp index a7fc3222a0c9..22b4026e9374 100644 --- a/UnitTests/SU2_CFD/periodic_limiters.cpp +++ b/UnitTests/SU2_CFD/periodic_limiters.cpp @@ -27,10 +27,25 @@ #include "catch.hpp" #include "../../SU2_CFD/include/solvers/CSolver.hpp" -TEST_CASE("Rotation of component bounds encloses every corner", "[PeriodicLimiter]") { - const su2double noAngles[3] = {0.0, 0.0, 0.0}, negativeAngle[3] = {0.0, -0.5, 0.0}; +TEST_CASE("Configured zero angles select translation and tiny angles select rotation", "[PeriodicLimiter]") { + const su2double noAngles[3] = {0.0, -0.0, 0.0}, negativeAngle[3] = {0.0, -0.5, 0.0}; CHECK_FALSE(GeometryToolbox::HasRotation(noAngles)); CHECK(GeometryToolbox::HasRotation(negativeAngle)); + const su2double deg2rad = PI_NUMBER / 180.0; + for (int iDim = 0; iDim < 3; ++iDim) { + su2double angles[3] = {0.0, 0.0, 0.0}; + angles[iDim] = su2double(0.0) * deg2rad; + CHECK_FALSE(GeometryToolbox::HasRotation(angles)); + angles[iDim] *= -1.0; + CHECK_FALSE(GeometryToolbox::HasRotation(angles)); + angles[iDim] = su2double(1.0e-18) * deg2rad; + CHECK(GeometryToolbox::HasRotation(angles)); + angles[iDim] *= -1.0; + CHECK(GeometryToolbox::HasRotation(angles)); + } +} + +TEST_CASE("Rotation of component bounds encloses every corner", "[PeriodicLimiter]") { auto checkCorners = [](auto& rotation, auto& lower, auto& upper) { constexpr int nDim = sizeof(lower) / sizeof(lower[0]); su2double expectedMin[nDim], expectedMax[nDim];