From bf973d4f7f386b71d3403b011ce7d790be0d23eb Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:33 +0200 Subject: [PATCH 1/6] Count the periodic volume and neighbors once with flow and heat solvers Both solvers accumulated the periodic volume and the number of neighbors on the same geometry, the heat solver now does it only without a flow solver. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- SU2_CFD/src/solvers/CHeatSolver.cpp | 16 ++-- .../chtPinArray_2d/periodic_weak_heat.cfg | 79 +++++++++++++++++++ TestCases/parallel_regression.py | 8 ++ 3 files changed, 97 insertions(+), 6 deletions(-) create mode 100644 TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index e23d25af07c2..302a80032b59 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -148,12 +148,16 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, const CSolver* fl ghostNodes = make_unique(Solution_Inf[0], maxMarkerVertices, nDim, nVar, config); } - /*--- Communicate and store volume and the number of neighbors for any dual CVs that lie on on periodic markers. ---*/ - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); + /*--- Communicate and store volume and the number of neighbors for any dual CVs that lie on on periodic markers. + * With a flow solver on the same geometry this was already done by the flow solver, and the values are + * accumulated, so it must not be done twice. ---*/ + if (!flow) { + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); + } } /*--- Store if implicit scheme is used. This has implications on the Residual and Jacobian handling for periodic * boundaries ---*/ diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg new file mode 100644 index 000000000000..6ab487b08981 --- /dev/null +++ b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg @@ -0,0 +1,79 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D pin array, periodic flow driven by a body force, with % +% the weakly coupled heat equation (flow and heat solver on % +% the same mesh, the heat solver does not change the flow) % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +SOLVER= INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO +% +INC_ENERGY_EQUATION= NO +WEAKLY_COUPLED_HEAT_EQUATION= YES + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_DENSITY_INIT= 1045.0 +INC_VELOCITY_INIT= ( 0.0001, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 338.0 +INC_NONDIM= DIMENSIONAL +% +SPECIFIC_HEAT_CP= 3540.0 +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.001385 +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +PRANDTL_LAM= 11.7 +% +BODY_FORCE= YES +BODY_FORCE_VECTOR= ( 0.5, 0.0, 0.0 ) + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_SYM= ( fluid_symmetry ) +MARKER_PERIODIC= ( fluid_inlet, fluid_outlet, 0.0,0.0,0.0, 0.0,0.0,0.0, 0.0111544,0.0,0.0 ) +MARKER_HEATFLUX= ( fluid_pin1_interface, 0.0, fluid_pin2_interface, 0.0, fluid_pin3_interface, 0.0 ) +% +MARKER_MONITORING= ( fluid_pin2_interface ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 1000 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1e-8 +LINEAR_SOLVER_ITER= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_HEAT= SCALAR_UPWIND +TIME_DISCRE_HEAT= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_PRESSURE +CONV_RESIDUAL_MINVAL= -16 +CONV_STARTITER= 10 +ITER= 2000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= fluid.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, DRAG ) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index e85b5b7be660..52254bcec9a9 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -780,6 +780,14 @@ def main(): inc_heatTransfer_BC.test_vals = [-8.904266, -7.745636, -8.064003, -1.079197, -1671.000000] test_list.append(inc_heatTransfer_BC) + # 2D pin array, periodic with a body force, flow and weakly coupled heat equation + inc_periodic_weak_heat = TestCase('inc_periodic_weak_heat') + inc_periodic_weak_heat.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d" + inc_periodic_weak_heat.cfg_file = "periodic_weak_heat.cfg" + inc_periodic_weak_heat.test_iter = 10 + inc_periodic_weak_heat.test_vals = [-4.712559, -5.747243, -5.751727, 873.702714] + test_list.append(inc_periodic_weak_heat) + ############################ ### Incompressible RANS ### ############################ From cd4b42b7b2fe7bcff940eb94b3ceb464a1fd05e6 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:33 +0200 Subject: [PATCH 2/6] Update the number of neighbors of halo points after the periodic count The periodic communication adds the neighbors of the periodic match to the owned points only, the halo copies kept the old count, which is used by the centered schemes. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 7 +++++++ TestCases/parallel_regression.py | 2 +- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index cb44757c6280..b68c9de9e5d2 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -260,6 +260,13 @@ void CFVMFlowSolverBase::CommunicateInitialState(CGeometry* geometry, cons InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); } + + /*--- The periodic communication updates the number of neighbors of the owned points only, update the halos. ---*/ + + if (config->GetnMarker_Periodic() > 0) { + geometry->InitiateComms(geometry, config, MPI_QUANTITIES::NEIGHBORS); + geometry->CompleteComms(geometry, config, MPI_QUANTITIES::NEIGHBORS); + } SetImplicitPeriodic(euler_implicit); if (MGLevel == MESH_0) SetRotatePeriodic(true); diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 52254bcec9a9..4565b4944eb7 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1172,7 +1172,7 @@ def main(): sbs_backward_step.cfg_dir = "backscatter/backward_step" sbs_backward_step.cfg_file = "backwardStep.cfg" sbs_backward_step.test_iter = 3 - sbs_backward_step.test_vals = [-6.352884, -3.465372, -5.507901, -3.906545, -9.506305, -6.365236, -6.331021, -6.331028] + sbs_backward_step.test_vals = [-6.352885, -3.465372, -5.507901, -3.906544, -9.506305, -6.365236, -6.331021, -6.331028] sbs_backward_step.unsteady = True sbs_backward_step.decompress = True sbs_backward_step.grid_file = "backward_step.su2" From 908a3c45b243c0c4dfc5013565f7dd1104a59c9a Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:33 +0200 Subject: [PATCH 3/6] Fix the pressure sensor of the MSW scheme at periodic points The division by the pressure was applied to the running maximum instead of the pressure difference of each edge. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- SU2_CFD/src/solvers/CHeatSolver.cpp | 2 +- SU2_CFD/src/solvers/CSolver.cpp | 2 +- TestCases/parallel_regression.py | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index 302a80032b59..2aacbf2cc4b3 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -152,7 +152,7 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, const CSolver* fl * With a flow solver on the same geometry this was already done by the flow solver, and the values are * accumulated, so it must not be done twice. ---*/ if (!flow) { - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) { + for (auto iPeriodic = 1u; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) { InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 56cee4b6ad52..fa87814c995d 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -703,7 +703,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, if ((!boundary_i || boundary_j) && geometry->nodes->GetDomain(iPoint)) { if (msw) { - Sensor_i = fmax(Sensor_i, fabs(Pressure_j - Pressure_i)) / fmin(Pressure_i, Pressure_j); + Sensor_i = fmax(Sensor_i, fabs(Pressure_j - Pressure_i) / fmin(Pressure_i, Pressure_j)); } else { Sensor_i += (Pressure_j - Pressure_i); Sensor_j += (Pressure_i + Pressure_j); diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 4565b4944eb7..3319ebf2f3e5 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -782,8 +782,8 @@ def main(): # 2D pin array, periodic with a body force, flow and weakly coupled heat equation inc_periodic_weak_heat = TestCase('inc_periodic_weak_heat') - inc_periodic_weak_heat.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d" - inc_periodic_weak_heat.cfg_file = "periodic_weak_heat.cfg" + inc_periodic_weak_heat.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d" + inc_periodic_weak_heat.cfg_file = "periodic_weak_heat.cfg" inc_periodic_weak_heat.test_iter = 10 inc_periodic_weak_heat.test_vals = [-4.712559, -5.747243, -5.751727, 873.702714] test_list.append(inc_periodic_weak_heat) From f289001b97d34c2428cc8623581116197bc61569 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:34:21 +0200 Subject: [PATCH 4/6] Complete intersecting periodic stencils and inherit coarse boundary flags --- Common/include/geometry/CGeometry.hpp | 27 +++ Common/src/geometry/CMultiGridGeometry.cpp | 9 +- .../computeGradientsLeastSquares.hpp | 2 + .../include/solvers/CFVMFlowSolverBase.hpp | 6 +- SU2_CFD/include/solvers/CSolver.hpp | 2 + SU2_CFD/src/solvers/CEulerSolver.cpp | 7 +- SU2_CFD/src/solvers/CSolver.cpp | 95 +++++++-- UnitTests/SU2_CFD/gradients.cpp | 193 ++++++++++++++++++ UnitTests/UnitQuadTestCase.hpp | 17 +- 9 files changed, 336 insertions(+), 22 deletions(-) diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 9d3e574b4328..1adcd124e4b7 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -62,6 +62,7 @@ extern "C" { #include "../CConfig.hpp" #include "../toolboxes/graph_toolbox.hpp" +#include "../toolboxes/geometry_toolbox.hpp" #include "../adt/CADTElemClass.hpp" using namespace std; @@ -517,6 +518,32 @@ class CGeometry { */ void AllocatePeriodicComms(unsigned short val_countPerPeriodicPoint); + /*! \brief Fraction of an edge stencil owned by this periodic copy. */ + inline passivedouble GetPeriodicEdgeWeight(unsigned long iPoint, unsigned long jPoint, const CConfig& config) const { + if (!nodes->GetPeriodicBoundary(iPoint) || !nodes->GetPeriodicBoundary(jPoint)) return 1.0; + passivedouble weight = 1.0; + for (auto iMarker = 0u; iMarker < nMarker; ++iMarker) { + if (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY || nodes->GetVertex(iPoint, iMarker) < 0 || + nodes->GetVertex(jPoint, iMarker) < 0) + continue; + const auto tag = config.GetMarker_All_TagBound(iMarker); + const auto donor = config.GetMarker_Periodic_Donor(tag); + if (donor < nMarker && nodes->GetVertex(iPoint, donor) >= 0 && nodes->GetVertex(jPoint, donor) >= 0) { + /*--- An edge along the rotation axis is shared by all N sectors, + * rather than by two copies of each of the two periodic faces. ---*/ + if (config.GetMarker_All_PerBound(iMarker) > config.GetnMarker_Periodic() / 2) continue; + const auto* angles = config.GetPeriodicRotAngles(tag); + su2double rotation[3][3]; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation); + const auto cosAngle = SU2_TYPE::GetValue(0.5 * (rotation[0][0] + rotation[1][1] + rotation[2][2] - 1)); + weight *= acos(std::max(-1.0, std::min(1.0, cosAngle))) / (2 * PI_NUMBER); + } else { + weight *= 0.5; + } + } + return weight; + } + /*! * \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/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index d3bca729bf5b..e02d613962ec 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1157,7 +1157,14 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); if (fine_grid->nodes->GetBoundary(iFinePoint)) { nodes->SetBoundary(iCoarsePoint, nMarker); - break; + nodes->SetPeriodicBoundary(iCoarsePoint, nodes->GetPeriodicBoundary(iCoarsePoint) || + fine_grid->nodes->GetPeriodicBoundary(iFinePoint)); + nodes->SetPhysicalBoundary(iCoarsePoint, nodes->GetPhysicalBoundary(iCoarsePoint) || + fine_grid->nodes->GetPhysicalBoundary(iFinePoint)); + nodes->SetSolidBoundary( + iCoarsePoint, nodes->GetSolidBoundary(iCoarsePoint) || fine_grid->nodes->GetSolidBoundary(iFinePoint)); + nodes->SetViscousBoundary( + iCoarsePoint, nodes->GetViscousBoundary(iCoarsePoint) || fine_grid->nodes->GetViscousBoundary(iFinePoint)); } } diff --git a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp index 296f433ff40f..a239f62989d6 100644 --- a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp +++ b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp @@ -253,6 +253,8 @@ void computeGradientsLeastSquares(CSolver* solver, if (weight > 0.0) { weight = 1.0 / weight; + if (config.GetnMarker_Periodic() > 2) + weight *= geometry.GetPeriodicEdgeWeight(iPoint, jPoint, config); for (size_t iDim = 0; iDim < nDim; ++iDim) for (size_t jDim = iDim; jDim < nDim; ++jDim) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index ae07bef84552..10e0bf538cbd 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -807,8 +807,10 @@ class CFVMFlowSolverBase : public CSolver { iPoint_UndLapl[iPoint] = fmax(iPoint_UndLapl[iPoint], fabs(sensVar_j - sensVar_i) / fmin(sensVar_j, sensVar_i)); } else { /*--- Jameson dissipation sensor, add variable difference and variable sum. ---*/ - iPoint_UndLapl[iPoint] += sensVar_j - sensVar_i; - jPoint_UndLapl[iPoint] += sensVar_j + sensVar_i; + const auto weight = config->GetnMarker_Periodic() > 2 ? + geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) : 1.0; + iPoint_UndLapl[iPoint] += weight * (sensVar_j - sensVar_i); + jPoint_UndLapl[iPoint] += weight * (sensVar_j + sensVar_i); } } diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 7a7b6c726021..feda8a8d6bb2 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -133,6 +133,8 @@ class CSolver { /*--- End variables that need to go. ---*/ + su2activevector periodicNeighborCount; /*!< \brief Partial periodic stencil count, including shared edges. */ + su2activevector iPoint_UndLapl; /*!< \brief Auxiliary variable for the undivided Laplacians. */ su2activevector jPoint_UndLapl; /*!< \brief Auxiliary variable for the undivided Laplacians. */ diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 32c265e55656..5f056d860473 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2404,13 +2404,16 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *co /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ if (boundary_i && !boundary_j) continue; + const auto weight = config->GetnMarker_Periodic() > 2 ? + geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) : 1.0; + /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ for (unsigned short iVar = 0; iVar < nVar; iVar++) - nodes->AddUnd_Lapl(iPoint, iVar, nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar)); + nodes->AddUnd_Lapl(iPoint, iVar, weight * (nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar))); su2double Pressure_j = nodes->GetPressure(jPoint); - nodes->AddUnd_Lapl(iPoint, nVar-1, Pressure_j-Pressure_i); + nodes->AddUnd_Lapl(iPoint, nVar-1, weight * (Pressure_j-Pressure_i)); } } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index fa87814c995d..e0c27ec7ac9d 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -207,7 +207,7 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, break; case PERIODIC_NEIGHBORS: COUNT_PER_POINT = 1; - MPI_TYPE = COMM_TYPE::UNSIGNED_SHORT; + MPI_TYPE = COMM_TYPE::DOUBLE; break; case PERIODIC_RESIDUAL: COUNT_PER_POINT = nVar + nVar*nVar + 1; @@ -352,6 +352,18 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, SU2_MPI::Error("The NEMO solvers do not support rotational periodicity yet.", CURRENT_FUNCTION); } + const bool intersectingPairs = config->GetnMarker_Periodic() > 2; + if (commType == PERIODIC_NEIGHBORS && intersectingPairs && val_periodic_index == 1) { + SU2_OMP_SAFE_GLOBAL_ACCESS(periodicNeighborCount.resize(geometry->GetnPoint());) + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) { + periodicNeighborCount(iPoint) = 0; + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) + periodicNeighborCount(iPoint) += geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config); + } + END_SU2_OMP_FOR + } + /*--- Local variables ---*/ bool boundary_i, boundary_j; @@ -413,7 +425,6 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, su2double *bufDSend = geometry->bufD_PeriodicSend; - unsigned short *bufSSend = geometry->bufS_PeriodicSend; /*--- Handle the different types of gradient and limiter. ---*/ @@ -501,6 +512,11 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, case PERIODIC_NEIGHBORS: + if (intersectingPairs) { + bufDSend[buf_offset] = periodicNeighborCount(iPoint); + break; + } + nNeighbor = 0; for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { @@ -508,13 +524,13 @@ 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 (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) nNeighbor++; } /*--- Store the number of neighbors in bufffer. ---*/ - bufSSend[buf_offset] = nNeighbor; + bufDSend[buf_offset] = nNeighbor; break; @@ -605,6 +621,16 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, case PERIODIC_LAPLACIAN: + if (intersectingPairs) { + for (auto iVar = 0u; iVar < nVar; ++iVar) + bufDSend[buf_offset + iVar] = base_nodes->GetUndivided_Laplacian(iPoint, iVar); + if (rotate_periodic) Rotate(zeros, &bufDSend[buf_offset + 1], &Und_Lapl[1]); + if (rotate_periodic) + for (auto iDim = 0u; iDim < nDim; ++iDim) + bufDSend[buf_offset + 1 + iDim] = Und_Lapl[1 + iDim]; + break; + } + /*--- For JST, the undivided Laplacian must be computed consistently by using the complete control volume info from both sides of the periodic face. ---*/ @@ -617,7 +643,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 (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) { /*--- Solution differences ---*/ @@ -671,6 +697,11 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, break; case PERIODIC_SENSOR: { + if (intersectingPairs) { + bufDSend[buf_offset] = iPoint_UndLapl(iPoint); + bufDSend[buf_offset + 1] = jPoint_UndLapl(iPoint); + break; + } const bool msw = config->GetKind_Upwind_Flow() == UPWIND::MSW; /*--- For the centered schemes, the sensor must be computed @@ -683,7 +714,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 (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) < 1.0) continue; /*--- Use density instead of pressure for incomp. flows. ---*/ @@ -781,6 +812,40 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Set a flag for unweighted or weighted least-squares. ---*/ + if (intersectingPairs) { + /*--- Rmatrix stores the upper triangle of the normal equations; + * entry (2,1) is a duplicate of (0,2) used by the LS factorization. ---*/ + su2double matrix[3][3] = {}, rotated[3][3] = {}; + for (auto iDim = 0u; iDim < nDim; ++iDim) + for (auto jDim = 0u; jDim < nDim; ++jDim) + matrix[iDim][jDim] = base_nodes->GetRmatrix(iPoint, min(iDim,jDim), max(iDim,jDim)); + for (auto iDim = 0u; iDim < nDim; ++iDim) + for (auto jDim = 0u; jDim < nDim; ++jDim) + for (auto kDim = 0u; kDim < nDim; ++kDim) + for (auto lDim = 0u; lDim < nDim; ++lDim) { + const auto qik = nDim == 2 ? rotMatrix2D[iDim][kDim] : rotMatrix3D[iDim][kDim]; + const auto qjl = nDim == 2 ? rotMatrix2D[jDim][lDim] : rotMatrix3D[jDim][lDim]; + rotated[iDim][jDim] += qik * matrix[kDim][lDim] * qjl; + } + for (auto iDim = 0u; iDim < nDim; ++iDim) + for (auto jDim = 0u; jDim < nDim; ++jDim) + bufDSend[buf_offset++] = iDim <= jDim ? rotated[iDim][jDim] : + (nDim == 3 && iDim == 2 && jDim == 1 ? rotated[0][2] : su2double(0)); + for (auto iVar = 0u; iVar < ICOUNT; ++iVar) + Rotate(zeros, gradient[iPoint][iVar], rotBlock[iVar]); + if (rotate_periodic) { + for (auto iDim = 0u; iDim < nDim; ++iDim) { + su2double velocity[3] = {}, rotatedVelocity[3] = {}; + for (auto jDim = 0u; jDim < nDim; ++jDim) velocity[jDim] = rotBlock[1+jDim][iDim]; + Rotate(zeros, velocity, rotatedVelocity); + for (auto jDim = 0u; jDim < nDim; ++jDim) rotBlock[1+jDim][iDim] = rotatedVelocity[jDim]; + } + } + for (auto iVar = 0u; iVar < ICOUNT; ++iVar) + for (auto iDim = 0u; iDim < nDim; ++iDim) bufDSend[buf_offset++] = rotBlock[iVar][iDim]; + break; + } + switch(commType) { case PERIODIC_SOL_ULS: case PERIODIC_SOL_ULS_R: @@ -826,7 +891,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 (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) { /*--- Get coordinates for the neighbor point. ---*/ @@ -1041,7 +1106,6 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, const su2double *bufDRecv = geometry->bufD_PeriodicRecv; - const unsigned short *bufSRecv = geometry->bufS_PeriodicRecv; /*--- Handle the different types of gradient and limiter. ---*/ @@ -1118,13 +1182,13 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, break; case PERIODIC_NEIGHBORS: - - /*--- Store the extra neighbors on the periodic face. ---*/ - - nNeighbor = (geometry->nodes->GetnNeighbor(iPoint) + - bufSRecv[buf_offset]); - geometry->nodes->SetnNeighbor(iPoint, nNeighbor); - + if (config->GetnMarker_Periodic() > 2) { + periodicNeighborCount(iPoint) += bufDRecv[buf_offset]; + geometry->nodes->SetnNeighbor(iPoint, SU2_TYPE::Int(periodicNeighborCount(iPoint) + 0.5)); + } else { + nNeighbor = geometry->nodes->GetnNeighbor(iPoint) + SU2_TYPE::Int(bufDRecv[buf_offset]); + geometry->nodes->SetnNeighbor(iPoint, nNeighbor); + } break; case PERIODIC_RESIDUAL: @@ -2280,6 +2344,7 @@ void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config) for (unsigned short iVar = 0; iVar < nVar; iVar++) { su2double delta = base_nodes->GetSolution(jPoint,iVar)-base_nodes->GetSolution(iPoint,iVar); + if (config->GetnMarker_Periodic() > 2) delta *= geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config); base_nodes->AddUnd_Lapl(iPoint, iVar, delta); } } diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index 72baf525ca68..671a8e1e262f 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -26,6 +26,8 @@ */ #include "catch.hpp" +#include "../UnitQuadTestCase.hpp" +#include "../../Common/include/geometry/CMultiGridGeometry.hpp" #include "../../Common/include/geometry/CPhysicalGeometry.hpp" #include "../../Common/include/containers/container_decorators.hpp" #include "../../SU2_CFD/include/solvers/CSolver.hpp" @@ -154,3 +156,194 @@ TEST_CASE("GG", "[Gradients]") { testGreenGauss(); } TEST_CASE("LS", "[Gradients]") { testLeastSquares(false); } TEST_CASE("WLS", "[Gradients]") { testLeastSquares(true); } + +TEST_CASE("Intersecting periodic stencils", "[Periodic][Gradients]") { + const auto nPairs = GENERATE(2u, 3u); + const bool diagonal = GENERATE(false, true); + const auto scheme = GENERATE(0u, 1u, 2u); + const bool msw = scheme == 1, incompressible = scheme == 2; + 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 == 2 ? "MARKER_CUSTOM= (z_plus,z_minus)\n" : ""); + std::string periodic = "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0, y_minus,y_plus, 0,0,0, 0,0,0, 0,1,0"; + if (nPairs == 3) periodic += ", z_minus,z_plus, 0,0,0, 0,0,0, 0,0,1"; + field.AddOption(periodic + ")"); + field.AddOption(msw ? "CONV_NUM_METHOD_FLOW= MSW" : "CONV_NUM_METHOD_FLOW= JST"); + field.AddOption("MUSCL_FLOW= NO\nNUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES"); + if (incompressible) { + field.SetOption("SOLVER= INC_NAVIER_STOKES"); + field.AddOption("FLUID_MODEL= CONSTANT_DENSITY"); + field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION"); + } + field.InitConfig(); + field.InitGeometry(true); + if (diagonal) { + /*--- A triangulated plane adds NE/SW edges. The SW donor at a periodic + * corner is reached through two successive pairs, rather than one. ---*/ + std::vector> neighbors(field.geometry->GetnPoint()); + for (auto i = 0ul; i < neighbors.size(); ++i) { + for (auto j : field.geometry->nodes->GetPoints(i)) neighbors[i].push_back(j); + const auto* x = field.geometry->nodes->GetCoord(i); + for (auto j = 0ul; j < neighbors.size(); ++j) { + const auto* y = field.geometry->nodes->GetCoord(j); + if (fabs(y[0] - x[0] - 0.25) < 1e-12 && fabs(y[1] - x[1] - 0.25) < 1e-12 && fabs(y[2] - x[2]) < 1e-12) { + neighbors[i].push_back(j); + neighbors[j].push_back(i); + } + } + } + field.geometry->nodes->SetPoints(neighbors); + delete field.geometry->edges; + field.geometry->SetEdges(); + for (auto i = 0ul; i < neighbors.size(); ++i) field.geometry->nodes->SetnNeighbor(i, neighbors[i].size()); + } + 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* solver = field.solver[FLOW_SOL]; + auto* nodes = solver->GetNodes(); + auto value = [](su2double x, su2double y, su2double z) { + return 10 + x * (1 - x) + 2 * y * (1 - y) + 3 * z * (1 - z); + }; + unsigned long target = field.geometry->GetnPoint(); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + nodes->SetPressure(i, value(x[0], x[1], x[2])); + for (auto v = 0u; v < solver->GetnVar(); ++v) + nodes->SetSolution( + i, v, + v == 0 ? value(x[0], x[1], x[2]) + : (v + 1 == solver->GetnVar() ? (incompressible ? 300.0 : value(x[0], x[1], x[2]) / 0.4) : 0)); + if (field.geometry->nodes->GetDomain(i) && x[0] == 0 && x[1] == 0 && x[2] == (nPairs == 3 ? 0 : 0.5)) target = i; + } + const auto z = nPairs == 3 ? 0.0 : 0.5; + const auto center = value(0, 0, z); + std::vector stencil = {value(0.25, 0, z), value(0.75, 0, z), + value(0, 0.25, z), value(0, 0.75, z), + value(0, 0, z + 0.25), value(0, 0, nPairs == 3 ? 0.75 : z - 0.25)}; + if (diagonal) { + stencil.push_back(value(0.25, 0.25, z)); + stencil.push_back(value(0.75, 0.75, z)); + } + su2double numerator = 0, denominator = 0, maximum = 0; + for (const auto neighbor : stencil) { + numerator += neighbor - center; + denominator += neighbor + center; + maximum = std::max(maximum, fabs(neighbor - center) / std::min(center, neighbor)); + } + field.config->SetGlobalParam(field.config->GetKind_Solver(), RUNTIME_FLOW_SYS); + { + SU2_OMP_PARALLEL { + solver->Preprocessing(field.geometry.get(), field.solver, field.config.get(), 0, 0, RUNTIME_FLOW_SYS, false); + } + END_SU2_OMP_PARALLEL + } + unsigned long localTarget = target < field.geometry->GetnPoint(), globalTarget = 0; + SU2_MPI::Allreduce(&localTarget, &globalTarget, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(globalTarget == 1); + INFO("pairs " << nPairs << ", diagonal " << diagonal << ", MSW " << msw << ", incompressible " << incompressible); + if (target < field.geometry->GetnPoint()) { + CHECK(field.geometry->nodes->GetnNeighbor(target) == stencil.size()); + CHECK(nodes->GetSensor(target) == + Approx(msw ? maximum : (incompressible ? 0.0 : fabs(numerator) / denominator)).margin(1e-12)); + if (!msw) CHECK(nodes->GetUndivided_Laplacian(target, 0) == Approx(numerator).margin(1e-12)); + } + /*--- Periodic scalar field: compare the assembled primitive LS normal + * equations with the complete stencil, including the composed diagonal. ---*/ + solver->SetRotatePeriodic(false); + auto& gradient = nodes->GetGradient_Primitive(); + auto& matrix = nodes->GetRmatrix(); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) nodes->SetPrimitive(i, v, value(x[0], x[1], x[2])); + } + SU2_OMP_PARALLEL { + computeGradientsLeastSquares(solver, MPI_QUANTITIES::PRIMITIVE_GRADIENT, PERIODIC_PRIM_LS, *field.geometry, + *field.config, true, nodes->GetPrimitive(), 0, solver->GetnPrimVarGrad(), -1, gradient, + matrix); + } + END_SU2_OMP_PARALLEL + if (target < field.geometry->GetnPoint()) { + CHECK(gradient(target, 0, 0) == Approx(0).margin(1e-12)); + CHECK(gradient(target, 0, 1) == Approx(0).margin(1e-12)); + CHECK(gradient(target, 0, 2) == Approx(nPairs == 3 ? 0 : 3 * (1 - 2 * z)).margin(1e-12)); + } +} + +TEST_CASE("Coarse periodic boundary flags", "[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_HEATFLUX= (y_minus,0,y_plus,0)\nMARKER_CUSTOM= (z_plus,z_minus)\n"); + field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); + field.AddOption("MGLEVEL= 1"); + field.InitConfig(); + field.InitGeometry(true); + CMultiGridGeometry coarse(field.geometry.get(), field.config.get(), 1); + coarse.SetPoint_Connectivity(field.geometry.get()); + coarse.SetVertex(field.geometry.get(), field.config.get()); + unsigned long periodic = 0, missing = 0; + for (auto marker = 0u; marker < coarse.GetnMarker(); ++marker) { + if (field.config->GetMarker_All_KindBC(marker) == HEAT_FLUX) { + for (auto v = 0ul; v < coarse.GetnVertex(marker); ++v) { + const auto i = coarse.vertex[marker][v]->GetNode(); + CHECK(coarse.nodes->GetPhysicalBoundary(i)); + CHECK(coarse.nodes->GetSolidBoundary(i)); + CHECK(coarse.nodes->GetViscousBoundary(i)); + } + } + if (field.config->GetMarker_All_KindBC(marker) != PERIODIC_BOUNDARY) continue; + for (auto vertex = 0ul; vertex < coarse.GetnVertex(marker); ++vertex) { + ++periodic; + if (!coarse.nodes->GetPeriodicBoundary(coarse.vertex[marker][vertex]->GetNode())) ++missing; + } + } + unsigned long local[] = {periodic, missing}, global[2] = {}; + SU2_MPI::Allreduce(local, global, 2, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(global[0] > 0); + CHECK(global[1] == 0); +} + +TEST_CASE("One-cell periodic stencil", "[Periodic][Gradients]") { + 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.SetOption("MESH_BOX_SIZE= 2,5,5"); + field.AddOption("NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES"); + field.InitConfig(); + field.InitGeometry(true); + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + field.InitSolver(); + auto* solver = field.solver[FLOW_SOL]; + auto* nodes = solver->GetNodes(); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) nodes->SetPrimitive(i, v, 3 + x[1] + x[2]); + } + solver->SetRotatePeriodic(false); + auto& gradient = nodes->GetGradient_Primitive(); + auto& matrix = nodes->GetRmatrix(); + SU2_OMP_PARALLEL { + computeGradientsLeastSquares(solver, MPI_QUANTITIES::PRIMITIVE_GRADIENT, PERIODIC_PRIM_LS, *field.geometry, + *field.config, true, nodes->GetPrimitive(), 0, solver->GetnPrimVarGrad(), -1, gradient, + matrix); + } + END_SU2_OMP_PARALLEL + unsigned long local = 0, global = 0; + for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + if (x[1] != 0.5 || x[2] != 0.5) continue; + ++local; + CHECK(field.geometry->nodes->GetnNeighbor(i) == 6); + CHECK(matrix(i, 0, 0) == Approx(2).margin(1e-12)); + CHECK(gradient(i, 0, 1) == Approx(1).margin(1e-12)); + CHECK(gradient(i, 0, 2) == Approx(1).margin(1e-12)); + } + SU2_MPI::Allreduce(&local, &global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(global == 2); +} diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp index d8da75dc0960..d943c1696ee5 100644 --- a/UnitTests/UnitQuadTestCase.hpp +++ b/UnitTests/UnitQuadTestCase.hpp @@ -61,6 +61,16 @@ struct UnitQuadTestCase { */ void AddOption(const std::string& optionLine) { config_options += optionLine + "\n"; } + /*! \brief Replace one existing base option without repeating its key. */ + void SetOption(const std::string& optionLine) { + const auto key = optionLine.substr(0, optionLine.find('=') + 1); + const auto start = config_options.find(key); + if (start == std::string::npos) + AddOption(optionLine); + else + config_options.replace(start, config_options.find('\n', start) - start, optionLine); + } + /*! * \brief Initialize the config structure */ @@ -83,10 +93,13 @@ struct UnitQuadTestCase { /*! * \brief Initialize the geometry */ - void InitGeometry() { + void InitGeometry(bool partition = false) { cout.rdbuf(nullptr); { - auto aux_geometry = std::unique_ptr(new CPhysicalGeometry(config.get(), 0, 1)); + const auto rank = partition ? SU2_MPI::GetRank() : 0; + const auto size = partition ? SU2_MPI::GetSize() : 1; + auto aux_geometry = std::unique_ptr(new CPhysicalGeometry(config.get(), rank, size)); + if (partition) aux_geometry->SetColorGrid_Parallel(config.get()); geometry = std::unique_ptr(new CPhysicalGeometry(aux_geometry.get(), config.get())); } geometry->SetSendReceive(config.get()); From c500a3ab3e2f67c09489351e36bb60241afbf4b6 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:53:45 +0200 Subject: [PATCH 5/6] Normalize assignment spacing reported by CodeFactor --- TestCases/parallel_regression.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 3319ebf2f3e5..b679fd1e7c13 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -781,7 +781,7 @@ def main(): test_list.append(inc_heatTransfer_BC) # 2D pin array, periodic with a body force, flow and weakly coupled heat equation - inc_periodic_weak_heat = TestCase('inc_periodic_weak_heat') + inc_periodic_weak_heat = TestCase('inc_periodic_weak_heat') inc_periodic_weak_heat.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d" inc_periodic_weak_heat.cfg_file = "periodic_weak_heat.cfg" inc_periodic_weak_heat.test_iter = 10 @@ -1172,7 +1172,7 @@ def main(): sbs_backward_step.cfg_dir = "backscatter/backward_step" sbs_backward_step.cfg_file = "backwardStep.cfg" sbs_backward_step.test_iter = 3 - sbs_backward_step.test_vals = [-6.352885, -3.465372, -5.507901, -3.906544, -9.506305, -6.365236, -6.331021, -6.331028] + sbs_backward_step.test_vals = [-6.352885, -3.465372, -5.507901, -3.906544, -9.506305, -6.365236, -6.331021, -6.331028] sbs_backward_step.unsteady = True sbs_backward_step.decompress = True sbs_backward_step.grid_file = "backward_step.su2" From afb78422e81e7492141bf5d0996aa5301ba46142 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 19:00:45 +0200 Subject: [PATCH 6/6] Avoid shadowing intersecting periodic communication indices --- SU2_CFD/src/solvers/CSolver.cpp | 40 +++++++++++++++++---------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index e0c27ec7ac9d..895da0e329e4 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -622,12 +622,12 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, case PERIODIC_LAPLACIAN: if (intersectingPairs) { - for (auto iVar = 0u; iVar < nVar; ++iVar) - bufDSend[buf_offset + iVar] = base_nodes->GetUndivided_Laplacian(iPoint, iVar); + for (auto iField = 0u; iField < nVar; ++iField) + bufDSend[buf_offset + iField] = base_nodes->GetUndivided_Laplacian(iPoint, iField); if (rotate_periodic) Rotate(zeros, &bufDSend[buf_offset + 1], &Und_Lapl[1]); if (rotate_periodic) - for (auto iDim = 0u; iDim < nDim; ++iDim) - bufDSend[buf_offset + 1 + iDim] = Und_Lapl[1 + iDim]; + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + bufDSend[buf_offset + 1 + iCoordinate] = Und_Lapl[1 + iCoordinate]; break; } @@ -816,33 +816,35 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Rmatrix stores the upper triangle of the normal equations; * entry (2,1) is a duplicate of (0,2) used by the LS factorization. ---*/ su2double matrix[3][3] = {}, rotated[3][3] = {}; - for (auto iDim = 0u; iDim < nDim; ++iDim) + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) for (auto jDim = 0u; jDim < nDim; ++jDim) - matrix[iDim][jDim] = base_nodes->GetRmatrix(iPoint, min(iDim,jDim), max(iDim,jDim)); - for (auto iDim = 0u; iDim < nDim; ++iDim) + matrix[iCoordinate][jDim] = + base_nodes->GetRmatrix(iPoint, min(iCoordinate,jDim), max(iCoordinate,jDim)); + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) for (auto jDim = 0u; jDim < nDim; ++jDim) for (auto kDim = 0u; kDim < nDim; ++kDim) for (auto lDim = 0u; lDim < nDim; ++lDim) { - const auto qik = nDim == 2 ? rotMatrix2D[iDim][kDim] : rotMatrix3D[iDim][kDim]; + const auto qik = nDim == 2 ? rotMatrix2D[iCoordinate][kDim] : rotMatrix3D[iCoordinate][kDim]; const auto qjl = nDim == 2 ? rotMatrix2D[jDim][lDim] : rotMatrix3D[jDim][lDim]; - rotated[iDim][jDim] += qik * matrix[kDim][lDim] * qjl; + rotated[iCoordinate][jDim] += qik * matrix[kDim][lDim] * qjl; } - for (auto iDim = 0u; iDim < nDim; ++iDim) + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) for (auto jDim = 0u; jDim < nDim; ++jDim) - bufDSend[buf_offset++] = iDim <= jDim ? rotated[iDim][jDim] : - (nDim == 3 && iDim == 2 && jDim == 1 ? rotated[0][2] : su2double(0)); - for (auto iVar = 0u; iVar < ICOUNT; ++iVar) - Rotate(zeros, gradient[iPoint][iVar], rotBlock[iVar]); + bufDSend[buf_offset++] = iCoordinate <= jDim ? rotated[iCoordinate][jDim] : + (nDim == 3 && iCoordinate == 2 && jDim == 1 ? rotated[0][2] : su2double(0)); + for (auto iField = 0u; iField < ICOUNT; ++iField) + Rotate(zeros, gradient[iPoint][iField], rotBlock[iField]); if (rotate_periodic) { - for (auto iDim = 0u; iDim < nDim; ++iDim) { + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) { su2double velocity[3] = {}, rotatedVelocity[3] = {}; - for (auto jDim = 0u; jDim < nDim; ++jDim) velocity[jDim] = rotBlock[1+jDim][iDim]; + for (auto jDim = 0u; jDim < nDim; ++jDim) velocity[jDim] = rotBlock[1+jDim][iCoordinate]; Rotate(zeros, velocity, rotatedVelocity); - for (auto jDim = 0u; jDim < nDim; ++jDim) rotBlock[1+jDim][iDim] = rotatedVelocity[jDim]; + for (auto jDim = 0u; jDim < nDim; ++jDim) rotBlock[1+jDim][iCoordinate] = rotatedVelocity[jDim]; } } - for (auto iVar = 0u; iVar < ICOUNT; ++iVar) - for (auto iDim = 0u; iDim < nDim; ++iDim) bufDSend[buf_offset++] = rotBlock[iVar][iDim]; + for (auto iField = 0u; iField < ICOUNT; ++iField) + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + bufDSend[buf_offset++] = rotBlock[iField][iCoordinate]; break; }