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 9a1489a77e18759141fb66ecda222457cf2bfde0 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:33 +0200 Subject: [PATCH 06/11] Treat periodic points on the rotation axis as points with N copies A point on the rotation axis is its own periodic match. It was handled like a pair of points: it received its own data rotated once in each direction, and in implicit runs its residual and Jacobian were then removed, so its solution never changed. Now it receives the data of every other copy of its control volume (the rotation applied 1 to N-1 times, N = 360 deg / angle), it keeps its equations, and after an implicit update its solution is averaged over the copies, which removes the velocity normal to the axis. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- Common/include/geometry/CGeometry.hpp | 10 +- Common/src/geometry/CGeometry.cpp | 100 +++++++++++++++--- SU2_CFD/src/solvers/CSolver.cpp | 92 ++++++++++++++-- .../navierstokes/periodic3D_axis/config.cfg | 69 ++++++++++++ TestCases/parallel_regression.py | 8 ++ 5 files changed, 254 insertions(+), 25 deletions(-) create mode 100644 TestCases/navierstokes/periodic3D_axis/config.cfg diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 9d3e574b4328..7d6a5f0b9bf1 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -342,8 +342,14 @@ class CGeometry { a particular vertex to be sent in periodic comms. */ *Local_Marker_PeriodicRecv{nullptr}; /*!< \brief Data structure holding the local index of the periodic marker for a particular vertex to be received in periodic comms. */ - su2double* bufD_PeriodicRecv{nullptr}; /*!< \brief Data structure for su2double periodic receive. */ - su2double* bufD_PeriodicSend{nullptr}; /*!< \brief Data structure for su2double periodic send. */ + vector Local_Copy_PeriodicSend; /*!< \brief For points on a rotation axis, which are their own periodic + match, the number of times the rotation is applied to the data + sent (one send per other copy of the control volume). Zero for + all other points. */ + vector Local_Copy_PeriodicRecv; /*!< \brief Same as Local_Copy_PeriodicSend, for the data received. */ + bool PeriodicAxisPoints{false}; /*!< \brief Whether this rank has periodic points on a rotation axis. */ + su2double* bufD_PeriodicRecv{nullptr}; /*!< \brief Data structure for su2double periodic receive. */ + su2double* bufD_PeriodicSend{nullptr}; /*!< \brief Data structure for su2double periodic send. */ unsigned short* bufS_PeriodicRecv{nullptr}; /*!< \brief Data structure for unsigned long periodic receive. */ unsigned short* bufS_PeriodicSend{nullptr}; /*!< \brief Data structure for unsigned long periodic send. */ SU2_MPI::Request* req_PeriodicSend{nullptr}; /*!< \brief Data structure for periodic send requests. */ diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 360da1aaaa1e..80d428eb6dea 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -1194,6 +1194,56 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { using PointMarkerSet = std::unordered_set; std::vector Points_Send_All(size, PointMarkerSet(0, pairHash)); + /*--- Points on a rotation axis are their own periodic match. The other copies of their control + volume are obtained by applying the rotation 1, 2, ... N-1 times (N = 360 deg / angle), instead + of once in each direction. Therefore, the vertex on the first marker of the pair sends N-1 times + to itself, and the vertex on the second marker does not send. ---*/ + + const unsigned long nMarker_All = config->GetnMarker_All(); + PeriodicAxisPoints = false; + + /*--- Number of copies of the domain around the rotation axis of a periodic marker (N). ---*/ + + auto nCopiesAroundAxis = [&](unsigned short jMarker) { + /*--- Tolerance, in radians, to accept the rotation angle as an integer fraction of 360 degrees. ---*/ + constexpr passivedouble angleTol = 1e-4; + + /*--- The trace of a rotation matrix is 1 + 2 cos(angle). ---*/ + const auto markerTag = config->GetMarker_All_TagBound(jMarker); + const su2double* angles = config->GetPeriodicRotAngles(markerTag); + su2double rotMatrix[MAXNDIM][MAXNDIM]; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotMatrix); + su2double cosAngle = 0.5 * (rotMatrix[0][0] + rotMatrix[1][1] + rotMatrix[2][2] - 1.0); + if (nDim == 2) cosAngle = cos(angles[2]); + cosAngle = min(cosAngle, 1.0); + cosAngle = max(cosAngle, -1.0); + const su2double angle = acos(cosAngle); + + const int nCopy = (angle > angleTol) ? SU2_TYPE::Int(2 * PI_NUMBER / angle + 0.5) : 0; + + if (fabs(nCopy * angle - 2 * PI_NUMBER) > angleTol) { + SU2_MPI::Error( + "Periodic points on the rotation axis require an angle that divides 360 degrees (marker " + markerTag + ").", + CURRENT_FUNCTION); + } + return nCopy; + }; + + /*--- Number of sends of a vertex (1 for the points that are not on a rotation axis), + and whether it is on a rotation axis. ---*/ + + auto nSendOfVertex = [&](unsigned short jMarker, unsigned long jVertex, bool& onAxis) -> unsigned long { + const auto* vert = geometry->vertex[jMarker][jVertex]; + onAxis = (static_cast(vert->GetDonorProcessor()) == rank) && + (static_cast(vert->GetDonorPoint()) == vert->GetNode()); + if (!onAxis) return 1; + if (config->GetMarker_All_PerBound(jMarker) > config->GetnMarker_Periodic() / 2) return 0; + + PeriodicAxisPoints = true; + return nCopiesAroundAxis(jMarker) - 1; + }; + bool onAxis = false; + /*--- Loop through all of our periodic markers and track our sends with each rank. ---*/ @@ -1213,8 +1263,11 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { iRank = static_cast(geometry->vertex[iMarker][iVertex]->GetDonorProcessor()); - /*--- Store the (point, marker) pair in the set for the destination rank. ---*/ - Points_Send_All[iRank].insert(std::make_pair(iPoint, static_cast(iMarker))); + /*--- Store the (point, marker) pair in the set for the destination rank, + once for each time the vertex sends (see above). ---*/ + const auto nCopy = nSendOfVertex(iMarker, iVertex, onAxis); + for (auto iCopy = 0ul; iCopy < nCopy; iCopy++) + Points_Send_All[iRank].insert(std::make_pair(iPoint, iMarker + iCopy * nMarker_All)); } } } @@ -1305,6 +1358,9 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { Local_Marker_PeriodicRecv = new unsigned long[nPoint_PeriodicRecv[nPeriodicRecv]]; for (iRecv = 0; iRecv < nPoint_PeriodicRecv[nPeriodicRecv]; iRecv++) Local_Marker_PeriodicRecv[iRecv] = 0; + Local_Copy_PeriodicSend.assign(nPoint_PeriodicSend[nPeriodicSend], 0); + Local_Copy_PeriodicRecv.assign(nPoint_PeriodicRecv[nPeriodicRecv], 0); + /*--- We allocate the buffers for communicating values in a later step once we know the maximum packet size that we need to communicate. This memory is deallocated and reallocated automatically in the case that @@ -1325,11 +1381,12 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { req_PeriodicRecv = new SU2_MPI::Request[nPeriodicRecv]; } - /*--- Allocate arrays for sending the periodic point index and marker - index to the recv rank so that it can store the local values. Therefore, - the recv rank can quickly loop through the buffers to unpack the data. ---*/ + /*--- Allocate arrays for sending the periodic point index, the marker index, + and the copy index (for points on a rotation axis) to the recv rank so that + it can store the local values. Therefore, the recv rank can quickly loop + through the buffers to unpack the data. ---*/ - unsigned short nPackets = 2; + const unsigned short nPackets = 3; auto* idSend = new unsigned long[nPoint_PeriodicSend[nPeriodicSend] * nPackets]; for (iSend = 0; iSend < nPoint_PeriodicSend[nPeriodicSend] * nPackets; iSend++) idSend[iSend] = 0; @@ -1366,17 +1423,24 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { if (iRank == destRank) { /*--- Check if we have already added this (point, marker) pair for this destination rank. Use the result if insert(), which is - a pair whose second element is success. ---*/ - const auto pointMarkerPair = std::make_pair(iPoint, static_cast(iMarker)); - const auto insertResult = Points_Send_All[destRank].insert(pointMarkerPair); - if (insertResult.second) { - Local_Point_PeriodicSend[ii] = iPoint; - Local_Marker_PeriodicSend[ii] = static_cast(iMarker); - jj = ii * nPackets; - idSend[jj] = geometry->vertex[iMarker][iVertex]->GetDonorPoint(); - jj++; - idSend[jj] = static_cast(iPeriodic); - ii++; + a pair whose second element is success. Points on a rotation + axis send more than once (see above). ---*/ + const auto nCopy = nSendOfVertex(iMarker, iVertex, onAxis); + for (auto iCopy = 0ul; iCopy < nCopy; iCopy++) { + const auto pointMarkerPair = std::make_pair(iPoint, iMarker + iCopy * nMarker_All); + const auto insertResult = Points_Send_All[destRank].insert(pointMarkerPair); + if (insertResult.second) { + Local_Point_PeriodicSend[ii] = iPoint; + Local_Marker_PeriodicSend[ii] = static_cast(iMarker); + Local_Copy_PeriodicSend[ii] = onAxis ? iCopy + 1 : 0; + jj = ii * nPackets; + idSend[jj] = geometry->vertex[iMarker][iVertex]->GetDonorPoint(); + jj++; + idSend[jj] = static_cast(iPeriodic); + jj++; + idSend[jj] = Local_Copy_PeriodicSend[ii]; + ii++; + } } } } @@ -1485,6 +1549,8 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { ii++; Local_Marker_PeriodicRecv[iRecv] = idRecv[ii]; ii++; + Local_Copy_PeriodicRecv[iRecv] = idRecv[ii]; + ii++; } delete[] idSend; diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index a8d90309eb3b..6ba1e8b1a698 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -342,6 +342,25 @@ namespace PeriodicCommHelpers { bool isRotation(const su2double* angles) { return angles[0] != 0.0 || angles[1] != 0.0 || angles[2] != 0.0; } + + /*--- Replaces a rotation matrix by the matrix that applies the same rotation nRot times. ---*/ + void applyRotationNTimes(unsigned long nRot, su2double rotMatrix[3][3]) { + su2double rotOnce[3][3], rotPrev[3][3]; + for (auto iDim = 0u; iDim < 3; iDim++) + for (auto jDim = 0u; jDim < 3; jDim++) + rotOnce[iDim][jDim] = rotMatrix[iDim][jDim]; + + for (auto iRot = 1ul; iRot < nRot; iRot++) { + for (auto iDim = 0u; iDim < 3; iDim++) + for (auto jDim = 0u; jDim < 3; jDim++) + rotPrev[iDim][jDim] = rotMatrix[iDim][jDim]; + + for (auto iDim = 0u; iDim < 3; iDim++) + for (auto jDim = 0u; jDim < 3; jDim++) + rotMatrix[iDim][jDim] = rotOnce[iDim][0]*rotPrev[0][jDim] + rotOnce[iDim][1]*rotPrev[1][jDim] + + rotOnce[iDim][2]*rotPrev[2][jDim]; + } + } } su2activematrix* CSolver::GetPeriodicProjections(const CConfig& config) { @@ -554,6 +573,32 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, GeometryToolbox::RotationMatrix(Theta, Phi, Psi, rotMatrix3D); } + /*--- Points on a rotation axis are their own periodic match, they send once for every + other copy of their control volume, i.e. with the rotation applied more than once. ---*/ + + const auto nRot = geometry->Local_Copy_PeriodicSend[msg_offset + iSend]; + + if (nRot > 1) { + if (nDim==2) { + const su2double PsiTotal = nRot*Psi; + GeometryToolbox::RotationMatrix(PsiTotal, rotMatrix2D); + } else { + PeriodicCommHelpers::applyRotationNTimes(nRot, rotMatrix3D); + } + } + + /*--- The neighbors on periodic faces are skipped where edges must be counted once, because + the periodic match also has these edges. From the second copy of a point on a rotation + axis onwards, the edges on the face of this marker (not along the axis) are new. ---*/ + + const auto donorMarker = (nRot > 1) ? config->GetMarker_Periodic_Donor(Marker_Tag) : 0; + + auto sharedEdge = [&](unsigned long jPoint) { + if (!geometry->nodes->GetPeriodicBoundary(jPoint)) return false; + return (nRot < 2) || (geometry->nodes->GetVertex(jPoint, iPeriodic) < 0) || + (geometry->nodes->GetVertex(jPoint, donorMarker) >= 0); + }; + /*--- Whether the vector components of the solution are rotated for this marker. ---*/ const bool rotation = rotate_periodic && PeriodicCommHelpers::isRotation(angles); @@ -587,7 +632,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, that we avoid double counting neighbors on both sides. If not, increment the count of neighbors for the donor. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) + if (!sharedEdge(jPoint)) nNeighbor++; } @@ -705,7 +750,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid periodic boundary points so that we do not duplicate edges on both sides of the periodic BC. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) { + if (!sharedEdge(jPoint)) { /*--- Solution differences ---*/ @@ -771,7 +816,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid halos and boundary points so that we don't duplicate edges on both sides of the periodic BC. ---*/ - if (geometry->nodes->GetPeriodicBoundary(jPoint)) continue; + if (sharedEdge(jPoint)) continue; /*--- Use density instead of pressure for incomp. flows. ---*/ @@ -914,7 +959,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid periodic boundary points so that we do not duplicate edges on both sides of the periodic BC. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) { + if (!sharedEdge(jPoint)) { /*--- Get coordinates for the neighbor point. ---*/ @@ -1204,7 +1249,10 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, nRecv = (geometry->nPoint_PeriodicRecv[jRecv+1] - geometry->nPoint_PeriodicRecv[jRecv]); - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + /*--- Points on a rotation axis receive more than once, in that case + the loop is not shared by the threads (one chunk). ---*/ + + SU2_OMP_FOR_STAT(geometry->PeriodicAxisPoints ? nRecv + 1 : size_t(OMP_MIN_SIZE)) for (iRecv = 0; iRecv < nRecv; iRecv++) { /*--- Get the local index for this communicated data. ---*/ @@ -1212,6 +1260,10 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, iPoint = geometry->Local_Point_PeriodicRecv[msg_offset + iRecv]; iPeriodic = geometry->Local_Marker_PeriodicRecv[msg_offset + iRecv]; + /*--- For points on a rotation axis, which copy of the control volume this is (0 otherwise). ---*/ + + const auto iCopy = geometry->Local_Copy_PeriodicRecv[msg_offset + iRecv]; + /*--- While all periodic face data was accumulated, we only store the values for the current pair of periodic faces. This is slightly inefficient when we have multiple pairs of periodic faces, but @@ -1283,6 +1335,18 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, Jacobian.AddBlock2Diag(iPoint, Jacobian_i); + /*--- A point on a rotation axis accumulates the residual and the diagonal block of all + the copies of its control volume, the blocks of its neighbors (which are the same in + all copies, up to the rotation) need the same factor, iCopy+1 after this copy. ---*/ + + if (iCopy > 0) { + const passivedouble factor = (iCopy + 1.0) / iCopy; + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + auto* block = Jacobian.GetBlock(iPoint, jPoint); + for (auto iEntry = 0; iEntry < nVar*nVar; iEntry++) block[iEntry] *= factor; + } + } + if (iPeriodic == val_periodic_index + nPeriodic/2) { for (iVar = 0; iVar < nVar; iVar++) { LinSysRes(iPoint, iVar) = 0.0; @@ -1317,6 +1381,22 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, } + /*--- A point on a rotation axis takes the average over all the copies of its control + volume (running average over the copies received), which removes the velocity + normal to the axis: a vector that is the same in all the copies is along the axis. ---*/ + + else if (iCopy > 0) { + + for (auto iVar = 0u; iVar < nVar; iVar++) { + const su2double average = base_nodes->GetSolution(iPoint, iVar) + + (bufDRecv[buf_offset] - base_nodes->GetSolution(iPoint, iVar)) / su2double(iCopy + 1); + base_nodes->SetSolution(iPoint, iVar, average); + base_nodes->SetSolution_Old(iPoint, iVar, average); + buf_offset++; + } + + } + break; case PERIODIC_LAPLACIAN: @@ -1324,7 +1404,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, /*--- Adjust the undivided Laplacian. The accumulation was with a subtraction before communicating, so now just add. ---*/ - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) base_nodes->AddUnd_Lapl(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; diff --git a/TestCases/navierstokes/periodic3D_axis/config.cfg b/TestCases/navierstokes/periodic3D_axis/config.cfg new file mode 100644 index 000000000000..fd1725ad14ac --- /dev/null +++ b/TestCases/navierstokes/periodic3D_axis/config.cfg @@ -0,0 +1,69 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Laminar flow in a pipe, 30 degree sector (one cell in the % +% circumferential direction) with nodes on the rotation axis % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +SOLVER= NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +% The initial flow is along the axis (z) +MACH_NUMBER= 0.2 +AOA= 90.0 +SIDESLIP_ANGLE= 0.0 +FREESTREAM_PRESSURE= 99000.0 +FREESTREAM_TEMPERATURE= 300.0 +REYNOLDS_NUMBER= 80 +REYNOLDS_LENGTH= 1.0 +REF_DIMENSIONALIZATION= DIMENSIONAL +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_PERIODIC= ( per1, per2, 0,0,0, 0,0,30, 0,0,0 ) +MARKER_HEATFLUX= ( wall, 0.0 ) +INLET_TYPE= TOTAL_CONDITIONS +MARKER_INLET= ( inlet, 300.0, 101000.0, 0.0, 0.0, 1.0 ) +MARKER_OUTLET= ( outlet, 100000.0 ) +% +MARKER_MONITORING= ( wall ) +MARKER_ANALYZE= ( outlet ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 10.0 +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1e-6 +LINEAR_SOLVER_ITER= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= NO + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_DENSITY +CONV_RESIDUAL_MINVAL= -12 +ITER= 5000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= wedge.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_DENSITY, RMS_MOMENTUM-Z, RMS_ENERGY, SURFACE_MASSFLOW ) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 87b0f125ff51..c3dc6a9074a8 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -445,6 +445,14 @@ def main(): periodic2d.test_vals = [-3.266248, -0.617784, -0.620252, 2.216476, -1.007761, -1035.000000] test_list.append(periodic2d) + # 3D rotational periodic pipe sector with nodes on the rotation axis + periodic3d_axis = TestCase('periodic3d_axis') + periodic3d_axis.cfg_dir = "navierstokes/periodic3D_axis" + periodic3d_axis.cfg_file = "config.cfg" + periodic3d_axis.test_iter = 100 + periodic3d_axis.test_vals = [-3.353165, 0.157501, 2.051458, -9.594400, -9.594400] + test_list.append(periodic3d_axis) + ########################## ### Compressible RANS ### ########################## From e0dfc6380fbd89faa666c35ee871e854e76b8892 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 08:28:39 +0200 Subject: [PATCH 07/11] Use the TestCases branch of this PR (revert before merging) The mesh of the new test periodic3d_axis is in the TestCases branch fix_periodic_axis. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- .github/workflows/regression.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index e1985a4911f2..8f08b2de4f6b 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -233,7 +233,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c fix_periodic_axis -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -282,7 +282,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests" + args: -b ${{github.ref}} -t develop -c fix_periodic_axis -s ${{matrix.testscript}} -a "--tapetests" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -330,7 +330,7 @@ jobs: PMIX_MCA_gds: hash with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan" + args: -b ${{github.ref}} -t develop -c fix_periodic_axis -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:260405-0054 with: @@ -375,7 +375,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan" + args: -b ${{github.ref}} -t develop -c fix_periodic_axis -s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: From 1925a8ef318b185ebf13bd09fb43222ff1f79f7b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:27:41 +0200 Subject: [PATCH 08/11] Complete periodic slip normals and refresh donor volumes after mesh updates --- Common/include/geometry/CGeometry.hpp | 6 + Common/src/geometry/CGeometry.cpp | 119 +++++++++++++++++- Common/src/geometry/CMultiGridGeometry.cpp | 2 + Common/src/geometry/CPhysicalGeometry.cpp | 2 + .../src/grid_movement/CVolumetricMovement.cpp | 2 + UnitTests/Common/geometry/CGeometry_test.cpp | 111 ++++++++++++++++ UnitTests/UnitQuadTestCase.hpp | 7 +- 7 files changed, 246 insertions(+), 3 deletions(-) diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 7d6a5f0b9bf1..f0cb724d7fe6 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -523,6 +523,12 @@ class CGeometry { */ void AllocatePeriodicComms(unsigned short val_countPerPeriodicPoint); + /*! \brief Sum partial scalar or vector geometry data over periodic copies (called by the master). */ + void SumPeriodicGeometry(const CConfig* config, su2activematrix& values, int vectorIndex); + + /*! \brief Refresh donor volumes after a mesh update. */ + void UpdatePeriodicVolumes(const CConfig* config); + /*! * \brief Routine to launch non-blocking recvs only for all periodic communication with neighboring partitions. * \note This routine is called by any class that has loaded data into the generic communication buffers. diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 80d428eb6dea..95055adc8c85 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -1555,6 +1555,77 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { delete[] idSend; delete[] idRecv; + ComputeModifiedSymmetryNormals(config); +} + +void CGeometry::SumPeriodicGeometry(const CConfig* config, su2activematrix& values, int vectorIndex) { + if (!nPeriodicSend && !nPeriodicRecv) return; + const auto count = values.cols(); + AllocatePeriodicComms(count); + for (auto iPair = 1u; iPair <= config->GetnMarker_Periodic() / 2; ++iPair) { + PostPeriodicRecvs(this, config, COMM_TYPE::DOUBLE, count); + for (auto iMessage = 0; iMessage < nPeriodicSend; ++iMessage) { + for (auto iSend = nPoint_PeriodicSend[iMessage]; iSend < nPoint_PeriodicSend[iMessage + 1]; ++iSend) { + const auto iPoint = Local_Point_PeriodicSend[iSend]; + auto* buffer = bufD_PeriodicSend + count * iSend; + for (auto iVar = 0u; iVar < count; ++iVar) buffer[iVar] = values(iPoint, iVar); + if (vectorIndex >= 0) { + const auto* angles = + config->GetPeriodicRotAngles(config->GetMarker_All_TagBound(Local_Marker_PeriodicSend[iSend])); + su2double rotation[3][3], vector[3] = {}; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation); + /*--- Axis points contribute the rotated normal of each additional copy. ---*/ + const auto nCopies = max(1ul, Local_Copy_PeriodicSend[iSend]); + for (auto iCopy = 0ul; iCopy < nCopies; ++iCopy) { + std::copy(buffer + vectorIndex, buffer + vectorIndex + nDim, vector); + for (auto iDim = 0u; iDim < nDim; ++iDim) { + buffer[vectorIndex + iDim] = 0; + for (auto jDim = 0u; jDim < nDim; ++jDim) + buffer[vectorIndex + iDim] += rotation[iDim][jDim] * vector[jDim]; + } + } + } + } +#ifdef HAVE_MPI + PostPeriodicSends(this, config, COMM_TYPE::DOUBLE, count, iMessage); +#else + const auto message = PeriodicRecv2Neighbor[rank]; + const auto start = count * nPoint_PeriodicSend[iMessage]; + const auto end = count * nPoint_PeriodicSend[iMessage + 1]; + std::copy(bufD_PeriodicSend + start, bufD_PeriodicSend + end, + bufD_PeriodicRecv + count * nPoint_PeriodicRecv[message]); +#endif + } + for (auto iMessage = 0; iMessage < nPeriodicRecv; ++iMessage) { + auto message = iMessage; +#ifdef HAVE_MPI + SU2_MPI::Status status; + int index; + SU2_MPI::Waitany(nPeriodicRecv, req_PeriodicRecv, &index, &status); + message = PeriodicRecv2Neighbor[status.MPI_SOURCE]; +#endif + for (auto iRecv = nPoint_PeriodicRecv[message]; iRecv < nPoint_PeriodicRecv[message + 1]; ++iRecv) { + const auto periodic = Local_Marker_PeriodicRecv[iRecv]; + if (periodic != iPair && periodic != iPair + config->GetnMarker_Periodic() / 2) continue; + const auto iPoint = Local_Point_PeriodicRecv[iRecv]; + for (auto iVar = 0u; iVar < count; ++iVar) values(iPoint, iVar) += bufD_PeriodicRecv[count * iRecv + iVar]; + } + } + SU2_MPI::Waitall(nPeriodicSend, req_PeriodicSend, MPI_STATUS_IGNORE); + } +} + +void CGeometry::UpdatePeriodicVolumes(const CConfig* config) { + if (!config->GetnMarker_Periodic() || (!nPeriodicSend && !nPeriodicRecv)) return; + AllocatePeriodicComms(1); + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + su2activematrix volumes(nPoint, 1); + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) volumes(iPoint, 0) = nodes->GetVolume(iPoint); + SumPeriodicGeometry(config, volumes, -1); + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) + nodes->SetPeriodicVolume(iPoint, volumes(iPoint, 0) - nodes->GetVolume(iPoint)); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS } void CGeometry::AllocatePeriodicComms(unsigned short countPerPeriodicPoint) { @@ -2800,6 +2871,7 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config) geometry_container[MESH_0]->SetControlVolume(config, UPDATE); geometry_container[MESH_0]->SetBoundControlVolume(config, UPDATE); geometry_container[MESH_0]->SetMaxLength(config); + geometry_container[MESH_0]->UpdatePeriodicVolumes(config); for (unsigned short iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { /*--- Update the control volume structures ---*/ @@ -2807,6 +2879,7 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config) geometry_container[iMesh]->SetControlVolume(geometry_container[iMesh - 1], UPDATE); geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh - 1], config, UPDATE); geometry_container[iMesh]->SetCoord(geometry_container[iMesh - 1]); + geometry_container[iMesh]->UpdatePeriodicVolumes(config); } /*--- Compute the global surface areas for all markers. ---*/ @@ -2943,6 +3016,46 @@ void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) { } } + std::vector>> periodicNormals(nMarker); + std::vector periodicPoints(nPoint, false); + if (nPeriodicRecv) + for (auto i = 0; i < nPoint_PeriodicRecv[nPeriodicRecv]; ++i) periodicPoints[Local_Point_PeriodicRecv[i]] = true; + + /*--- A slip wall meeting a periodic face needs the normal of the complete + * wall patch. Keep markers separate to preserve distinct corner constraints. ---*/ + if (config->GetnMarker_Periodic() && (nPeriodicSend || nPeriodicRecv)) { + for (auto iCfgMarker = 0u; iCfgMarker < config->GetnMarker_CfgFile(); ++iCfgMarker) { + const auto tag = config->GetMarker_CfgFile_TagBound(iCfgMarker); + const auto kind = config->GetMarker_CfgFile_KindBC(tag); + if (kind != SYMMETRY_PLANE && kind != EULER_WALL) continue; + su2activematrix normals(nPoint, nDim); + normals = su2double(0); + for (const auto iMarker : symMarkers) { + if (config->GetMarker_All_TagBound(iMarker) != tag) continue; + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; ++iVertex) { + const auto iPoint = vertex[iMarker][iVertex]->GetNode(); + for (auto iDim = 0u; iDim < nDim; ++iDim) normals(iPoint, iDim) = vertex[iMarker][iVertex]->GetNormal(iDim); + } + } + SumPeriodicGeometry(config, normals, 0); + for (const auto iMarker : symMarkers) { + if (config->GetMarker_All_TagBound(iMarker) != tag) continue; + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; ++iVertex) { + const auto iPoint = vertex[iMarker][iVertex]->GetNode(); + if (!periodicPoints[iPoint]) continue; + const auto area = GeometryToolbox::Norm(nDim, normals[iPoint]); + if (area <= MIN_AREA) continue; + auto& normal = symmetryNormals[iMarker][iVertex]; + auto& areaNormal = periodicNormals[iMarker][iVertex]; + for (auto iDim = 0u; iDim < nDim; ++iDim) { + areaNormal[iDim] = normals(iPoint, iDim); + normal[iDim] = areaNormal[iDim] / area; + } + } + } + } + } + /*--- Merge the normals of curved markers to be independent of how surfaces are divided. ---*/ std::unordered_map> mergedNormals; @@ -2959,7 +3072,11 @@ void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) { if (count < 2) continue; std::array normal = {}; - vertex[iMarker][iVertex]->GetNormal(normal.data()); + const auto corrected = periodicNormals[iMarker].find(iVertex); + if (corrected != periodicNormals[iMarker].end()) + normal = corrected->second; + else + vertex[iMarker][iVertex]->GetNormal(normal.data()); auto result = mergedNormals.emplace(iPoint, normal); const auto inserted = result.second; diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index d3bca729bf5b..5b3a2e58977a 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1333,6 +1333,8 @@ void CMultiGridGeometry::SetBoundControlVolume(const CGeometry* fine_grid, const } END_SU2_OMP_FOR + /*--- Allocate on the whole team before the master computes periodic slip normals. ---*/ + if (nPeriodicSend || nPeriodicRecv) AllocatePeriodicComms(nDim); SU2_OMP_SAFE_GLOBAL_ACCESS(ComputeModifiedSymmetryNormals(config);) } diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 6928c2d922b5..18af5b585ca1 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -6984,6 +6984,8 @@ void CPhysicalGeometry::SetBoundControlVolume(const CConfig* config, unsigned sh } END_SU2_OMP_FOR + /*--- Allocate on the whole team before the master computes periodic slip normals. ---*/ + if (nPeriodicSend || nPeriodicRecv) AllocatePeriodicComms(nDim); SU2_OMP_SAFE_GLOBAL_ACCESS(ComputeModifiedSymmetryNormals(config);) } diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp index ca6001b7f9f9..4e7677ea040c 100644 --- a/Common/src/grid_movement/CVolumetricMovement.cpp +++ b/Common/src/grid_movement/CVolumetricMovement.cpp @@ -45,6 +45,7 @@ dual mesh control volumes in the domain and on the boundaries. ---*/ geometry->SetControlVolume(config, UPDATE); geometry->SetBoundControlVolume(config, UPDATE); geometry->SetMaxLength(config); + geometry->UpdatePeriodicVolumes(config); } void CVolumetricMovement::UpdateMultiGrid(CGeometry** geometry, CConfig* config) { @@ -58,6 +59,7 @@ including computing the grid velocities on the coarser levels. ---*/ geometry[iMGlevel]->SetControlVolume(geometry[iMGfine], UPDATE); geometry[iMGlevel]->SetBoundControlVolume(geometry[iMGfine], config, UPDATE); geometry[iMGlevel]->SetCoord(geometry[iMGfine]); + geometry[iMGlevel]->UpdatePeriodicVolumes(config); if (config->GetGrid_Movement()) geometry[iMGlevel]->SetRestricted_GridVelocity(geometry[iMGfine]); } } diff --git a/UnitTests/Common/geometry/CGeometry_test.cpp b/UnitTests/Common/geometry/CGeometry_test.cpp index 22556b16055d..e2ff2695a208 100644 --- a/UnitTests/Common/geometry/CGeometry_test.cpp +++ b/UnitTests/Common/geometry/CGeometry_test.cpp @@ -26,7 +26,9 @@ */ #include "catch.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../UnitQuadTestCase.hpp" +#include "../../../Common/include/geometry/CMultiGridGeometry.hpp" std::unique_ptr TestCase; @@ -140,3 +142,112 @@ 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 slip-wall normal", "[Periodic]") { + const bool multigrid = GENERATE(false, true); + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, "MARKER_EULER= (y_minus,y_plus)\nMARKER_CUSTOM= (z_minus,z_plus)\n"); + field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 30,0,0, 0,0,0)"); + if (multigrid) field.AddOption("MGLEVEL= 1"); + field.InitConfig(); + field.InitGeometry(true); + /*--- Bend the BOX into an annular sector about x: y is radius, x is azimuth. ---*/ + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + const su2double angle = x[0] * PI_NUMBER / 6, radius = 1 + x[1], axial = x[2]; + const su2double coordinate[] = {axial, -radius * sin(angle), radius * cos(angle)}; + field.geometry->nodes->SetCoord(i, coordinate); + } + field.geometry->SetControlVolume(field.config.get(), UPDATE); + field.geometry->SetBoundControlVolume(field.config.get(), UPDATE); + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + std::unique_ptr coarse; + CGeometry* geometry = field.geometry.get(); + if (multigrid) { + coarse.reset(new CMultiGridGeometry(geometry, field.config.get(), 1)); + coarse->SetPoint_Connectivity(geometry); + coarse->SetEdges(); + coarse->SetVertex(geometry, field.config.get()); + coarse->SetControlVolume(geometry, ALLOCATE); + coarse->SetBoundControlVolume(geometry, field.config.get(), ALLOCATE); + coarse->SetCoord(geometry); + coarse->SetMGLevel(1); + coarse->MatchPeriodic(field.config.get(), 1); + coarse->PreprocessPeriodicComms(coarse.get(), field.config.get()); + geometry = coarse.get(); + } + su2double error = 0; + unsigned long checked = 0; + for (auto marker = 0u; marker < geometry->GetnMarker(); ++marker) { + if (field.config->GetMarker_All_TagBound(marker) != "y_plus") continue; + for (auto v = 0ul; v < geometry->GetnVertex(marker); ++v) { + const auto i = geometry->vertex[marker][v]->GetNode(); + const auto* x = geometry->nodes->GetCoord(i); + if (fabs(x[1]) > 1e-12 || fabs(x[2] - 2) > 1e-12 || !geometry->nodes->GetDomain(i)) continue; + const auto it = geometry->symmetryNormals[marker].find(v); + const auto* normal = + it == geometry->symmetryNormals[marker].end() ? geometry->vertex[marker][v]->GetNormal() : it->second.data(); + error = std::max(error, fabs(normal[1]) / GeometryToolbox::Norm(3, normal)); + ++checked; + } + } + unsigned long total = 0; + SU2_MPI::Allreduce(&checked, &total, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(total > 0); + CHECK(error < 1e-12); +} + +TEST_CASE("Periodic volume refresh after mesh deformation", "[Periodic]") { + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n"); + field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); + field.AddOption("DEFORM_MESH= YES"); + field.InitConfig(); + field.InitGeometry(true); + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + field.InitSolver(); + std::vector> original(field.geometry->GetnPoint()); + for (auto i = 0ul; i < original.size(); ++i) + for (auto d = 0u; d < 3; ++d) original[i][d] = field.geometry->nodes->GetCoord(i, d); + CGeometry* meshes[] = {field.geometry.get()}; + for (const auto amplitude : {0.05, -0.02, 0.0}) { + for (auto i = 0ul; i < original.size(); ++i) { + const auto& x = original[i]; + field.geometry->nodes->SetCoord(i, 0, x[0] + amplitude * sin(2 * PI_NUMBER * x[0]) * sin(PI_NUMBER * x[2])); + } + SU2_OMP_PARALLEL { CGeometry::UpdateGeometry(meshes, field.config.get()); } + END_SU2_OMP_PARALLEL + su2double expected = 0, stored = 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] - 1) < 1e-12) expected = field.geometry->nodes->GetVolume(i); + if (fabs(x[0]) < 1e-12) stored = field.geometry->nodes->GetPeriodicVolume(i); + } + su2double totals[2] = {stored, expected}, global[2] = {}; + SU2_MPI::Allreduce(totals, global, 2, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(global[1] > 0); + CHECK(global[0] == Approx(global[1]).margin(1e-12)); + auto* solver = field.solver[FLOW_SOL]; + auto* nodes = solver->GetNodes(); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) + for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) + nodes->SetPrimitive(i, v, 2 + field.geometry->nodes->GetCoord(i, 1)); + SU2_OMP_PARALLEL { solver->SetPrimitive_Gradient_GG(field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + su2double error = 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 && fabs(x[0] - 1) > 1e-12) continue; + error = std::max(error, fabs(nodes->GetGradient_Primitive()(i, 0, 1) - 1)); + } + CHECK(error < 1e-12); + } +} diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp index d8da75dc0960..b6e8bfdd6bc3 100644 --- a/UnitTests/UnitQuadTestCase.hpp +++ b/UnitTests/UnitQuadTestCase.hpp @@ -83,10 +83,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()); From a99faf7cb180a4e8a43c251af069b9b052b33250 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:36:11 +0200 Subject: [PATCH 09/11] Complete periodic implicit products, solves, and reverse coupling --- Common/include/linear_algebra/CSysMatrix.hpp | 11 ++ Common/src/linear_algebra/CSysMatrix.cpp | 45 ++++- Common/src/linear_algebra/CSysSolve.cpp | 15 ++ .../include/solvers/CFVMFlowSolverBase.hpp | 3 +- SU2_CFD/include/solvers/CScalarSolver.inl | 3 +- SU2_CFD/include/solvers/CSolver.hpp | 10 +- SU2_CFD/src/solvers/CSolver.cpp | 26 ++- .../edge_residual_blocks_tests.cpp | 169 ++++++++++++++++++ .../linear_algebra/periodic_ad_tests.cpp | 103 +++++++++++ UnitTests/UnitQuadTestCase.hpp | 10 ++ UnitTests/meson.build | 1 + 11 files changed, 385 insertions(+), 11 deletions(-) create mode 100644 UnitTests/Common/linear_algebra/periodic_ad_tests.cpp diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 70dd943690ec..dc621086b732 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -257,6 +257,13 @@ class CSysMatrix { private: friend struct CSysMatrixComms; + int periodicVectorIndex{-2}; /*!< \brief -2 disables periodic projection; -1 denotes scalar variables. */ + mutable CSysVector projectedInput; + mutable su2activematrix periodicBuffer; + + void ProjectPeriodic(const CSysVector& input, CSysVector& output, CGeometry* geometry, + const CConfig* config) const; + const int rank; /*!< \brief MPI Rank. */ const int size; /*!< \brief MPI Size. */ @@ -700,6 +707,10 @@ class CSysMatrix { void ComputeLU_SGSPreconditionerBackward(CSysVector& prod) const; public: + /*! \brief Keep all partial periodic rows and couple their matching unknowns in the product. */ + void SetPeriodicProjection(int vectorIndex) { periodicVectorIndex = vectorIndex; } + bool HasPeriodicProjection() const { return periodicVectorIndex >= -1; } + /*! * \brief Constructor of the class. */ diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 503e2ad4893c..74ec213541db 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -944,11 +944,41 @@ void CSysMatrix::DeleteValsRowi(unsigned long block_i, unsigned long } } +template +void CSysMatrix::ProjectPeriodic(const CSysVector& input, CSysVector& output, + CGeometry* geometry, const CConfig* config) const { + geometry->AllocatePeriodicComms(nVar + 1); + SU2_OMP_SAFE_GLOBAL_ACCESS(periodicBuffer.resize(nPoint, nVar + 1);) + SU2_OMP_FOR_STAT(omp_heavy_size) + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + for (auto iVar = 0u; iVar < nVar; ++iVar) periodicBuffer(iPoint, iVar) = input(iPoint, iVar); + periodicBuffer(iPoint, nVar) = 1; + } + END_SU2_OMP_FOR + SU2_OMP_SAFE_GLOBAL_ACCESS(geometry->SumPeriodicGeometry(config, periodicBuffer, periodicVectorIndex);) + SU2_OMP_FOR_STAT(omp_heavy_size) + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) + for (auto iVar = 0u; iVar < nVar; ++iVar) + output(iPoint, iVar) = ActiveAssign(periodicBuffer(iPoint, iVar) / periodicBuffer(iPoint, nVar)); + END_SU2_OMP_FOR + CSysMatrixComms::Initiate(output, geometry, config); + CSysMatrixComms::Complete(output, geometry, config); +} + template void CSysMatrix::MatrixVectorProduct(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const { SU2_ZONE_SCOPED + const bool periodic = HasPeriodicProjection() && config->GetnMarker_Periodic(); + if (periodic) { + if (projectedInput.GetLocSize() != nPoint * nVar) { + SU2_OMP_SAFE_GLOBAL_ACCESS(projectedInput.Initialize(nPoint, nPointDomain, nVar, ScalarType(0));) + } + ProjectPeriodic(vec, projectedInput, geometry, config); + } + const auto& input = periodic ? projectedInput : vec; + if (useCuda) { #ifdef SU2_ENABLE_CUDA_KERNELS if constexpr (su2_gpu_capable_v) { @@ -983,17 +1013,28 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector& v if (quantized_mode) { SU2_OMP_FOR_DYN(omp_heavy_size) for (auto row_i = 0ul; row_i < nPointDomain; row_i++) { - QuantizedRowProduct(vec, row_i, &prod[row_i * nVar]); + QuantizedRowProduct(input, row_i, &prod[row_i * nVar]); } END_SU2_OMP_FOR } else { SU2_OMP_FOR_DYN(omp_heavy_size) for (auto row_i = 0ul; row_i < nPointDomain; row_i++) { - RowProduct(vec, row_i, &prod[row_i * nVar]); + RowProduct(input, row_i, &prod[row_i * nVar]); } END_SU2_OMP_FOR } + if (periodic) { + /*--- K = P A P + I-P. P is an orthogonal average of periodic copies, + * so transposing A also gives the transpose of the complete operator. + * I-P enforces consistency without introducing missing-neighbour rows. ---*/ + ProjectPeriodic(prod, prod, geometry, config); + SU2_OMP_FOR_STAT(omp_heavy_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) + for (auto iVar = 0u; iVar < nVar; ++iVar) prod(iPoint, iVar) += vec(iPoint, iVar) - projectedInput(iPoint, iVar); + END_SU2_OMP_FOR + } + /*--- MPI Parallelization. ---*/ CSysMatrixComms::Initiate(prod, geometry, config); diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 8ec1b19b23b3..e1ea5fa855c9 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -45,6 +45,17 @@ SU2_RESTORE_WARNING #include namespace { +/*--- Both forward and reverse solves must use the periodic host product. ---*/ +void checkPeriodicSolver(bool projection, unsigned short kindSolver, const CConfig* config) { + if (!projection || !config->GetnMarker_Periodic()) return; + if (config->GetCUDA()) SU2_MPI::Error("Implicit periodic coupling requires ENABLE_CUDA= NO.", CURRENT_FUNCTION); + if (kindSolver == PASTIX_LU || kindSolver == PASTIX_LDLT) + SU2_MPI::Error( + "Direct PaStiX solves do not assemble periodic constraints. Use FGMRES with " + "LINEAR_SOLVER_PREC= PASTIX_LU instead.", + CURRENT_FUNCTION); +} + /*! * \brief Epsilon used in CSysSolve depending on datatype to * decide if the linear system is already solved. @@ -1451,6 +1462,8 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con } } + checkPeriodicSolver(Jacobian.HasPeriodicProjection(), KindSolver, config); + const bool nested = SetupInnerSolver(KindSolver, config); /*--- Stop the recording for the linear solver ---*/ @@ -1675,6 +1688,8 @@ unsigned long CSysSolve::Solve_b(CSysMatrix& Jacobian, c } } + checkPeriodicSolver(Jacobian.HasPeriodicProjection(), KindSolver, config); + const bool nested = SetupInnerSolver(KindSolver, config); /*--- Set up preconditioner and matrix-vector product ---*/ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index ae07bef84552..1d7ad96e9c59 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -979,7 +979,8 @@ class CFVMFlowSolverBase : public CSolver { if (nodes->GetDelta_Time(iPoint) != 0.0) { - su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); + su2double Vol = geometry->nodes->GetVolume(iPoint); + if (!Jacobian.HasPeriodicProjection()) Vol += geometry->nodes->GetPeriodicVolume(iPoint); su2double Delta = Vol / nodes->GetDelta_Time(iPoint); diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index af03659fe7f5..2e6b875c57b4 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -498,7 +498,8 @@ void CScalarSolver::PrepareImplicitIteration(CGeometry* geometry, const su2double dt = nodes->GetDelta_Time(iPoint); if (dt != 0.0) { - su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); + su2double Vol = geometry->nodes->GetVolume(iPoint); + if (!Jacobian.HasPeriodicProjection()) Vol += geometry->nodes->GetPeriodicVolume(iPoint); Jacobian.AddVal2Diag(iPoint, Vol / dt); } else { Jacobian.SetVal2Diag(iPoint, 1.0); diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 88e09d166dd1..6473f5d893d5 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -4224,13 +4224,19 @@ class CSolver { * \brief Routine that sets the flag controlling implicit treatment for periodic BCs. * \param[in] val_implicit_periodic - Flag controlling implicit treatment for periodic BCs. */ - inline void SetImplicitPeriodic(bool val_implicit_periodic) { implicit_periodic = val_implicit_periodic; } + inline void SetImplicitPeriodic(bool val_implicit_periodic) { + implicit_periodic = val_implicit_periodic; + Jacobian.SetPeriodicProjection(val_implicit_periodic ? (rotate_periodic ? 1 : -1) : -2); + } /*! * \brief Routine that sets the flag controlling solution rotation for periodic BCs. * \param[in] val_implicit_periodic - Flag controlling solution rotation for periodic BCs. */ - inline void SetRotatePeriodic(bool val_rotate_periodic) { rotate_periodic = val_rotate_periodic; } + inline void SetRotatePeriodic(bool val_rotate_periodic) { + rotate_periodic = val_rotate_periodic; + if (implicit_periodic) Jacobian.SetPeriodicProjection(val_rotate_periodic ? 1 : -1); + } /*! * \brief Storage for the limiters with rotational periodicity: the min and max, over the edges of the periodic diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 6ba1e8b1a698..af01bd8be729 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -210,7 +210,7 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, MPI_TYPE = COMM_TYPE::UNSIGNED_SHORT; break; case PERIODIC_RESIDUAL: - COUNT_PER_POINT = nVar + nVar*nVar + 1; + COUNT_PER_POINT = Jacobian.HasPeriodicProjection() ? nVar + 1 : nVar + nVar*nVar + 1; MPI_TYPE = COMM_TYPE::DOUBLE; break; case PERIODIC_IMPLICIT: @@ -667,7 +667,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, contributions to the Jacobian block diagonal, i.e., the impact of the point upon itself, J_ii. ---*/ - if (implicit_periodic) { + if (implicit_periodic && !Jacobian.HasPeriodicProjection()) { const auto block = Jacobian.GetBlockView(iPoint, iPoint); @@ -1198,7 +1198,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, su2double Time_Step, Volume; su2double **Jacobian_i = nullptr; - if ((commType == PERIODIC_RESIDUAL) && implicit_periodic) { + if ((commType == PERIODIC_RESIDUAL) && implicit_periodic && !Jacobian.HasPeriodicProjection()) { Jacobian_i = new su2double* [nVar]; for (iVar = 0; iVar < nVar; iVar++) Jacobian_i[iVar] = new su2double [nVar]; @@ -1305,7 +1305,13 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, /*--- Add contributions to total residual. ---*/ - LinSysRes.AddBlock(iPoint, &bufDRecv[buf_offset]); + if (Jacobian.HasPeriodicProjection()) { + const auto copies = iCopy > 0 ? iCopy : 1ul; + for (auto iVar = 0u; iVar < nVar; ++iVar) + LinSysRes(iPoint, iVar) = (copies * LinSysRes(iPoint, iVar) + bufDRecv[buf_offset + iVar]) / (copies + 1); + } else { + LinSysRes.AddBlock(iPoint, &bufDRecv[buf_offset]); + } buf_offset += nVar; /*--- Check the computed time step against the donor @@ -1324,7 +1330,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, the passive face such that it does not participate in the linear solve. ---*/ - if (implicit_periodic) { + if (implicit_periodic && !Jacobian.HasPeriodicProjection()) { for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { @@ -1367,6 +1373,16 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, we are updating the solution at the passive nodes using the new solution from the master. ---*/ + if (implicit_periodic && Jacobian.HasPeriodicProjection()) { + const auto copies = iCopy > 0 ? iCopy : 1ul; + for (auto iVar = 0u; iVar < nVar; ++iVar) { + const auto value = (copies * base_nodes->GetSolution(iPoint, iVar) + bufDRecv[buf_offset + iVar]) / (copies + 1); + base_nodes->SetSolution(iPoint, iVar, value); + base_nodes->SetSolution_Old(iPoint, iVar, value); + } + break; + } + if ((implicit_periodic) && (iPeriodic == val_periodic_index + nPeriodic/2)) { diff --git a/UnitTests/Common/linear_algebra/edge_residual_blocks_tests.cpp b/UnitTests/Common/linear_algebra/edge_residual_blocks_tests.cpp index 12cd5e04220b..453b01314fe4 100644 --- a/UnitTests/Common/linear_algebra/edge_residual_blocks_tests.cpp +++ b/UnitTests/Common/linear_algebra/edge_residual_blocks_tests.cpp @@ -26,6 +26,7 @@ */ #include "catch.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../UnitQuadTestCase.hpp" /*--- A block whose entries are all different, so a mixed-up index reads back wrong. ---*/ @@ -165,3 +166,171 @@ TEST_CASE("SetBlocks and SetOffDiagBlocks with quantized off-diagonal storage", CheckBlock(matrix, jPoint, iPoint, nVar, jac_ji_new, quantTol); } } + +TEST_CASE("Complete periodic implicit operator and transpose", "[Periodic][LinearAlgebra]") { + const auto kind = GENERATE(0u, 1u, 2u, 3u); + const bool rotation = kind == 1; + const auto nPairs = kind > 1 ? kind : 1u; + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, + nPairs == 1 ? "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n" + : (nPairs == 2 ? "MARKER_CUSTOM= (z_plus,z_minus)\n" : "")); + std::string periodic = rotation ? "MARKER_PERIODIC= (x_minus,x_plus, 0,0.5,0.5, 90,0,0, 1,0,0" + : "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0"; + if (nPairs > 1) periodic += ", y_minus,y_plus, 0,0,0, 0,0,0, 0,1,0"; + if (nPairs > 2) periodic += ", z_minus,z_plus, 0,0,0, 0,0,0, 0,0,1"; + field.AddOption(periodic + ")"); + /*--- Avoid asking a float Krylov solver to converge below roundoff. ---*/ + field.AddOption("LINEAR_SOLVER= BCGSTAB\nLINEAR_SOLVER_PREC= JACOBI\nLINEAR_SOLVER_ITER= 150"); + field.AddOption(sizeof(su2mixedfloat) == sizeof(float) ? "LINEAR_SOLVER_ERROR= 1e-6" : "LINEAR_SOLVER_ERROR= 1e-10"); + field.InitConfig(); + field.InitGeometry(true); + for (auto pair = 1u; pair <= nPairs; ++pair) field.geometry->MatchPeriodic(field.config.get(), pair); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + field.InitSolver(); + auto& matrix = field.solver[FLOW_SOL]->Jacobian; + const auto nVar = field.solver[FLOW_SOL]->GetnVar(); + const auto nPoint = field.geometry->GetnPoint(); + const auto nDomain = field.geometry->GetnPointDomain(); + unsigned long globalPoints = 0; + SU2_MPI::Allreduce(&nDomain, &globalPoints, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(nDomain > 0); + const auto size = globalPoints * nVar; + std::vector partialA(size * size, 0), partialP(size * size, 0), A(size * size), P(size * size); + matrix.SetValZero(); + for (auto i = 0ul; i < nPoint; ++i) { + const auto globalI = field.geometry->nodes->GetGlobalIndex(i); + std::vector columns = {i}; + for (auto j : field.geometry->nodes->GetPoints(i)) columns.push_back(j); + for (auto j : columns) { + const auto globalJ = field.geometry->nodes->GetGlobalIndex(j); + auto* block = matrix.GetBlock(i, j); + REQUIRE(block != nullptr); + for (auto a = 0u; a < nVar; ++a) + for (auto b = 0u; b < nVar; ++b) { + /*--- Nonsymmetric blocks expose mistakes in the reverse operator. ---*/ + const auto value = + (i == j && a == b ? 3.0 : 0.0) + 0.001 * (1 + a + 2 * b) + (i == j ? 0 : 0.002 * (1 + globalI)); + block[a * nVar + b] = value; + if (i < nDomain) + partialA[(globalI * nVar + a) * size + globalJ * nVar + b] = SU2_TYPE::GetValue(block[a * nVar + b]); + } + } + if (i < nDomain) + for (auto a = 0u; a < nVar; ++a) partialP[(globalI * nVar + a) * size + globalI * nVar + a] = 1; + } + SU2_MPI::Allreduce(partialA.data(), A.data(), A.size(), MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(partialP.data(), P.data(), P.size(), MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + for (auto pair = 1u; pair <= nPairs; ++pair) { + std::fill(partialP.begin(), partialP.end(), 0); + for (auto i = 0ul; i < nDomain; ++i) { + const auto globalI = field.geometry->nodes->GetGlobalIndex(i); + for (auto a = 0u; a < nVar; ++a) + std::copy_n(P.data() + (globalI * nVar + a) * size, size, partialP.data() + (globalI * nVar + a) * size); + } + for (auto marker = 0u; marker < field.geometry->GetnMarker(); ++marker) { + if (field.config->GetMarker_All_KindBC(marker) != PERIODIC_BOUNDARY) continue; + const auto index = static_cast(field.config->GetMarker_All_PerBound(marker)); + if (index != pair && index != pair + nPairs) continue; + const auto* angles = field.config->GetPeriodicRotAngles(field.config->GetMarker_All_TagBound(marker)); + su2double q[3][3]; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], q); + for (auto vertex = 0ul; vertex < field.geometry->GetnVertex(marker); ++vertex) { + const auto* point = field.geometry->vertex[marker][vertex]; + const auto i = point->GetNode(); + if (i >= nDomain) continue; + const auto globalI = field.geometry->nodes->GetGlobalIndex(i); + const auto donor = point->GetDonorGlobalIndex(); + for (auto a = 0u; a < nVar; ++a) + for (auto j = 0ul; j < size; ++j) { + auto value = P[(globalI * nVar + a) * size + j]; + if (a >= 1 && a <= 3) { + for (auto b = 1u; b <= 3; ++b) value += q[b - 1][a - 1] * P[(donor * nVar + b) * size + j]; + } else { + value += P[(donor * nVar + a) * size + j]; + } + partialP[(globalI * nVar + a) * size + j] = 0.5 * value; + } + } + } + SU2_MPI::Allreduce(partialP.data(), P.data(), P.size(), MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + } + auto product = [&](const std::vector& mat, const std::vector& x, bool transpose = false) { + std::vector result(size, 0); + for (auto i = 0ul; i < size; ++i) + for (auto j = 0ul; j < size; ++j) result[i] += (transpose ? mat[j * size + i] : mat[i * size + j]) * x[j]; + return result; + }; + std::vector x(size), y(size); + for (auto i = 0ul; i < size; ++i) { + x[i] = sin(0.1 * i); + y[i] = cos(0.07 * i); + } + const auto px = product(P, x), py = product(P, y); + auto reference = product(P, product(A, px)); + auto transposeReference = product(P, product(A, py, true)); + for (auto i = 0ul; i < size; ++i) { + reference[i] += x[i] - px[i]; + transposeReference[i] += y[i] - py[i]; + } + CSysVector input(nPoint, nDomain, nVar), output(nPoint, nDomain, nVar), + transposeOutput(nPoint, nDomain, nVar); + for (auto i = 0ul; i < nPoint; ++i) + for (auto a = 0u; a < nVar; ++a) input(i, a) = x[field.geometry->nodes->GetGlobalIndex(i) * nVar + a]; + SU2_OMP_PARALLEL { matrix.MatrixVectorProduct(input, output, field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + su2double error = 0; + for (auto i = 0ul; i < nDomain; ++i) + for (auto a = 0u; a < nVar; ++a) + error = std::max(error, fabs(SU2_TYPE::GetValue(output(i, a)) - + reference[field.geometry->nodes->GetGlobalIndex(i) * nVar + a])); + CHECK(error < 1e-5); + CSysSolve system; + CSysVector rhs(nPoint, nDomain, nVar), solution(nPoint, nDomain, nVar); + solution = su2double(0); + for (auto i = 0ul; i < nPoint; ++i) + for (auto a = 0u; a < nVar; ++a) rhs(i, a) = reference[field.geometry->nodes->GetGlobalIndex(i) * nVar + a]; + SU2_OMP_PARALLEL { system.Solve(matrix, rhs, solution, field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + error = 0; + for (auto i = 0ul; i < nDomain; ++i) + for (auto a = 0u; a < nVar; ++a) + error = std::max( + error, fabs(SU2_TYPE::GetValue(solution(i, a)) - x[field.geometry->nodes->GetGlobalIndex(i) * nVar + a])); + CHECK(error < 1e-5); + SU2_OMP_PARALLEL { matrix.TransposeInPlace(); } + END_SU2_OMP_PARALLEL + for (auto i = 0ul; i < nPoint; ++i) + for (auto a = 0u; a < nVar; ++a) input(i, a) = y[field.geometry->nodes->GetGlobalIndex(i) * nVar + a]; + SU2_OMP_PARALLEL { matrix.MatrixVectorProduct(input, transposeOutput, field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + error = 0; + su2double dotForward = 0, dotTranspose = 0; + for (auto i = 0ul; i < nDomain; ++i) + for (auto a = 0u; a < nVar; ++a) { + const auto global = field.geometry->nodes->GetGlobalIndex(i) * nVar + a; + error = std::max(error, fabs(SU2_TYPE::GetValue(transposeOutput(i, a)) - transposeReference[global])); + dotForward += y[global] * SU2_TYPE::GetValue(output(i, a)); + dotTranspose += x[global] * SU2_TYPE::GetValue(transposeOutput(i, a)); + } + CHECK(error < 1e-5); + su2double dots[] = {dotForward, dotTranspose}, globalDots[2] = {}; + SU2_MPI::Allreduce(dots, globalDots, 2, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + CHECK(globalDots[0] == Approx(globalDots[1]).margin(1e-4)); + SU2_OMP_PARALLEL { matrix.TransposeInPlace(); } + END_SU2_OMP_PARALLEL + solution = su2double(0); + for (auto i = 0ul; i < nPoint; ++i) + for (auto a = 0u; a < nVar; ++a) + rhs(i, a) = transposeReference[field.geometry->nodes->GetGlobalIndex(i) * nVar + a]; + SU2_OMP_PARALLEL { system.Solve_b(matrix, rhs, solution, field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + error = 0; + for (auto i = 0ul; i < nDomain; ++i) + for (auto a = 0u; a < nVar; ++a) + error = std::max( + error, fabs(SU2_TYPE::GetValue(solution(i, a)) - y[field.geometry->nodes->GetGlobalIndex(i) * nVar + a])); + CHECK(error < 1e-5); +} diff --git a/UnitTests/Common/linear_algebra/periodic_ad_tests.cpp b/UnitTests/Common/linear_algebra/periodic_ad_tests.cpp new file mode 100644 index 000000000000..79582cd69853 --- /dev/null +++ b/UnitTests/Common/linear_algebra/periodic_ad_tests.cpp @@ -0,0 +1,103 @@ +/*! + * \file periodic_ad_tests.cpp + * \brief Finite-difference check of the periodic linear-solve external adjoint. + * \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 "../../UnitQuadTestCase.hpp" + +TEST_CASE("Periodic linear solve external adjoint", "[Periodic][AD tests]") { + const bool rotation = GENERATE(false, true); + AD::Reset(); + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n"); + field.SetOption("SOLVER= EULER"); + field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION"); + field.AddOption("MATH_PROBLEM= DISCRETE_ADJOINT"); + field.AddOption(rotation ? "MARKER_PERIODIC= (x_minus,x_plus, 0,0.5,0.5, 90,0,0, 1,0,0)" + : "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); + field.AddOption("LINEAR_SOLVER= FGMRES\nLINEAR_SOLVER_PREC= JACOBI\nLINEAR_SOLVER_ERROR= 1e-12"); + field.AddOption("LINEAR_SOLVER_ITER= 200\nDISCADJ_LIN_SOLVER= FGMRES\nDISCADJ_LIN_PREC= JACOBI"); + field.InitConfig(); + field.InitGeometry(true); + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + const auto nPoint = field.geometry->GetnPoint(), nDomain = field.geometry->GetnPointDomain(); + REQUIRE(nDomain > 0); + constexpr unsigned short nVar = 5; + CSysMatrix matrix; + matrix.Initialize(nPoint, nDomain, nVar, nVar, true, field.geometry.get(), field.config.get()); + matrix.SetPeriodicProjection(1); + matrix.SetValZero(); + for (auto i = 0ul; i < nPoint; ++i) { + std::vector columns = {i}; + for (auto j : field.geometry->nodes->GetPoints(i)) columns.push_back(j); + for (auto j : columns) { + auto* block = matrix.GetBlock(i, j); + for (auto a = 0u; a < nVar; ++a) + for (auto b = 0u; b < nVar; ++b) + block[a * nVar + b] = (i == j && a == b ? 3.0 : 0.0) + 0.001 * (1 + a + 2 * b) + + (i == j ? 0 : 0.002 * (1 + field.geometry->nodes->GetGlobalIndex(i))); + } + } + CSysSolve system; + CSysVector rhs(nPoint, nDomain, nVar), solution(nPoint, nDomain, nVar); + auto solve = [&](const su2double& parameter) { + rhs = su2double(0); + solution = su2double(0); + for (auto i = 0ul; i < nDomain; ++i) + for (auto a = 0u; a < nVar; ++a) { + const auto global = field.geometry->nodes->GetGlobalIndex(i) * nVar + a; + rhs(i, a) = cos(0.03 * global) + parameter * sin(0.1 * global); + } + system.Solve(matrix, rhs, solution, field.geometry.get(), field.config.get()); + su2double objective = 0; + for (auto i = 0ul; i < nDomain; ++i) + for (auto a = 0u; a < nVar; ++a) + objective += cos(0.07 * (field.geometry->nodes->GetGlobalIndex(i) * nVar + a)) * solution(i, a); + return objective; + }; + su2double parameter = 0.2; + AD::StartRecording(); + AD::RegisterInput(parameter); + auto objective = solve(parameter); + AD::RegisterOutput(objective); + AD::StopRecording(); + SU2_TYPE::SetDerivative(objective, 1.0); + AD::ComputeAdjoint(); + su2double localDerivative = SU2_TYPE::GetDerivative(parameter), derivative = 0; + SU2_MPI::Allreduce(&localDerivative, &derivative, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + /*--- The external solve transposes its matrix during the reverse sweep. ---*/ + matrix.TransposeInPlace(); + AD::Reset(); + constexpr double step = 1e-5; + const auto plus = SU2_TYPE::GetValue(solve(su2double(0.2 + step))); + const auto minus = SU2_TYPE::GetValue(solve(su2double(0.2 - step))); + su2double localDifference = (plus - minus) / (2 * step), difference = 0; + SU2_MPI::Allreduce(&localDifference, &difference, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + CHECK(SU2_TYPE::GetValue(derivative) == Approx(SU2_TYPE::GetValue(difference)).epsilon(1e-6).margin(1e-6)); + AD::Reset(); +} diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp index b6e8bfdd6bc3..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 */ diff --git a/UnitTests/meson.build b/UnitTests/meson.build index ba3c63afce92..40b0ce4d7bc6 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -28,6 +28,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', # Reverse-mode (algorithmic differentiation) tests: su2_cfd_tests_ad = files(['Common/simple_ad_test.cpp', + 'Common/linear_algebra/periodic_ad_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests_AD.cpp', 'Common/toolboxes/multilayer_perceptron/MLP_Jacobian_tests.cpp']) From 9f9b20cf61332f555cd1ecadd205931082e1a0c4 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:38:37 +0200 Subject: [PATCH 10/11] Use precision-aware Krylov checks in the periodic operator testcase --- .../linear_algebra/edge_residual_blocks_tests.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/UnitTests/Common/linear_algebra/edge_residual_blocks_tests.cpp b/UnitTests/Common/linear_algebra/edge_residual_blocks_tests.cpp index 453b01314fe4..bc4706b76da1 100644 --- a/UnitTests/Common/linear_algebra/edge_residual_blocks_tests.cpp +++ b/UnitTests/Common/linear_algebra/edge_residual_blocks_tests.cpp @@ -183,8 +183,8 @@ TEST_CASE("Complete periodic implicit operator and transpose", "[Periodic][Linea if (nPairs > 2) periodic += ", z_minus,z_plus, 0,0,0, 0,0,0, 0,0,1"; field.AddOption(periodic + ")"); /*--- Avoid asking a float Krylov solver to converge below roundoff. ---*/ - field.AddOption("LINEAR_SOLVER= BCGSTAB\nLINEAR_SOLVER_PREC= JACOBI\nLINEAR_SOLVER_ITER= 150"); - field.AddOption(sizeof(su2mixedfloat) == sizeof(float) ? "LINEAR_SOLVER_ERROR= 1e-6" : "LINEAR_SOLVER_ERROR= 1e-10"); + field.AddOption("LINEAR_SOLVER_PREC= JACOBI\nLINEAR_SOLVER_ITER= 150"); + field.AddOption(sizeof(su2mixedfloat) == sizeof(float) ? "LINEAR_SOLVER_ERROR= 1e-7" : "LINEAR_SOLVER_ERROR= 1e-10"); field.InitConfig(); field.InitGeometry(true); for (auto pair = 1u; pair <= nPairs; ++pair) field.geometry->MatchPeriodic(field.config.get(), pair); @@ -287,6 +287,9 @@ TEST_CASE("Complete periodic implicit operator and transpose", "[Periodic][Linea error = std::max(error, fabs(SU2_TYPE::GetValue(output(i, a)) - reference[field.geometry->nodes->GetGlobalIndex(i) * nVar + a])); CHECK(error < 1e-5); + /*--- Krylov solution accuracy depends on the arithmetic used by its basis. ---*/ + const auto solveTolerance = + sizeof(su2mixedfloat) == sizeof(float) ? sqrt(std::numeric_limits::epsilon()) : 1e-5; CSysSolve system; CSysVector rhs(nPoint, nDomain, nVar), solution(nPoint, nDomain, nVar); solution = su2double(0); @@ -299,7 +302,7 @@ TEST_CASE("Complete periodic implicit operator and transpose", "[Periodic][Linea for (auto a = 0u; a < nVar; ++a) error = std::max( error, fabs(SU2_TYPE::GetValue(solution(i, a)) - x[field.geometry->nodes->GetGlobalIndex(i) * nVar + a])); - CHECK(error < 1e-5); + CHECK(error < solveTolerance); SU2_OMP_PARALLEL { matrix.TransposeInPlace(); } END_SU2_OMP_PARALLEL for (auto i = 0ul; i < nPoint; ++i) @@ -332,5 +335,5 @@ TEST_CASE("Complete periodic implicit operator and transpose", "[Periodic][Linea for (auto a = 0u; a < nVar; ++a) error = std::max( error, fabs(SU2_TYPE::GetValue(solution(i, a)) - y[field.geometry->nodes->GetGlobalIndex(i) * nVar + a])); - CHECK(error < 1e-5); + CHECK(error < solveTolerance); } From 4db1d1c3f6babb639f3e66b51f8d259c4b5aaf60 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:53:46 +0200 Subject: [PATCH 11/11] Normalize assignment spacing reported by CodeFactor --- TestCases/hybrid_regression.py | 2 +- TestCases/parallel_regression.py | 4 ++-- TestCases/serial_regression.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index cc24dd1f2601..baf87d0ac02f 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -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.634571, -8.895558, -9.348754] + 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 c3dc6a9074a8..44d7ce7acf88 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -446,7 +446,7 @@ def main(): test_list.append(periodic2d) # 3D rotational periodic pipe sector with nodes on the rotation axis - periodic3d_axis = TestCase('periodic3d_axis') + periodic3d_axis = TestCase('periodic3d_axis') periodic3d_axis.cfg_dir = "navierstokes/periodic3D_axis" periodic3d_axis.cfg_file = "config.cfg" periodic3d_axis.test_iter = 100 @@ -1309,7 +1309,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.634558, -8.895554, -9.348754] + 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 e7c0261da35f..8a59db7f9b90 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -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.634558, -8.895554, -9.348754] + 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)