From bf973d4f7f386b71d3403b011ce7d790be0d23eb Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:33 +0200 Subject: [PATCH 01/16] 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 02/16] 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 03/16] 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 04/16] 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 05/16] 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 06/16] 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; } From e5a3370d4ffd85acac36ee936e9e5d7800fe49d3 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 18:35:52 +0200 Subject: [PATCH 07/16] Cache periodic edge weights during geometry preprocessing --- Common/include/geometry/CGeometry.hpp | 32 +++++------------ Common/src/geometry/CGeometry.cpp | 32 +++++++++++++++++ .../computeGradientsLeastSquares.hpp | 7 ++-- .../include/solvers/CFVMFlowSolverBase.hpp | 7 ++-- SU2_CFD/src/solvers/CEulerSolver.cpp | 6 ++-- SU2_CFD/src/solvers/CSolver.cpp | 35 ++++++++++++------- UnitTests/SU2_CFD/gradients.cpp | 33 +++++++++++++++++ 7 files changed, 108 insertions(+), 44 deletions(-) diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 8d0af501cdd9..7ab090845020 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -62,7 +62,6 @@ extern "C" { #include "../CConfig.hpp" #include "../toolboxes/graph_toolbox.hpp" -#include "../toolboxes/geometry_toolbox.hpp" #include "../adt/CADTElemClass.hpp" using namespace std; @@ -320,6 +319,8 @@ class CGeometry { PassiveRequest* reqP_P2PSend{nullptr}; /*!< \brief Data structure for point-to-point send requests. */ PassiveRequest* reqP_P2PRecv{nullptr}; /*!< \brief Data structure for point-to-point recv requests. */ + vector periodicEdgeWeight; /*!< Fraction of each edge stencil owned by this periodic copy. */ + /*--- Data structures for periodic communications. ---*/ int maxCountPerPeriodicPoint{0}; /*!< \brief Maximum number of pieces of data sent per vertex in periodic comms. */ @@ -520,30 +521,13 @@ class CGeometry { */ void AllocatePeriodicComms(unsigned short val_countPerPeriodicPoint); + /*! \brief Rebuild after changing connectivity, periodic membership or transforms; not needed for coordinate motion. + */ + void SetPeriodicEdgeWeights(const CConfig& config); + /*! \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; + inline passivedouble GetPeriodicEdgeWeight(unsigned long iEdge) const { + return periodicEdgeWeight.empty() ? 1.0 : periodicEdgeWeight[iEdge]; } /*! diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index acafc2b7c025..400cb6e99beb 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -1155,7 +1155,38 @@ void CGeometry::MatchPeriodic(const CConfig* config, unsigned short val_periodic delete[] Buffer_Recv_Marker; } +void CGeometry::SetPeriodicEdgeWeights(const CConfig& config) { + periodicEdgeWeight.assign(nEdge, 1.0); + for (auto iEdge = 0ul; iEdge < nEdge; ++iEdge) { + const auto iPoint = edges->GetNode(iEdge, 0), jPoint = edges->GetNode(iEdge, 1); + if (!nodes->GetPeriodicBoundary(iPoint) || !nodes->GetPeriodicBoundary(jPoint)) continue; + 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; + } + } + periodicEdgeWeight[iEdge] = weight; + } +} + void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { + geometry->SetPeriodicEdgeWeights(*config); + /*--- We start with the send and receive lists already available in the form of stored periodic point-donor pairs. We will loop through these markers and establish the neighboring ranks and number of @@ -1648,6 +1679,7 @@ void CGeometry::PostPeriodicSends(CGeometry* geometry, const CConfig* config, CO } void CGeometry::SetEdges() { + periodicEdgeWeight.clear(); nEdge = 0; for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { for (auto iNode = 0u; iNode < nodes->GetnPoint(iPoint); iNode++) { diff --git a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp index a239f62989d6..eb259db87132 100644 --- a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp +++ b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp @@ -231,8 +231,9 @@ void computeGradientsLeastSquares(CSolver* solver, Rmatrix(iPoint, iDim, jDim) = 0.0; - for (auto jPoint : nodes->GetPoints(iPoint)) - { + for (auto iNeigh = 0u; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = nodes->GetEdge(iPoint, iNeigh); const auto coord_j = geometry.nodes->GetCoord(jPoint); AD::SetPreaccIn(coord_j, nDim); @@ -254,7 +255,7 @@ void computeGradientsLeastSquares(CSolver* solver, { weight = 1.0 / weight; if (config.GetnMarker_Periodic() > 2) - weight *= geometry.GetPeriodicEdgeWeight(iPoint, jPoint, config); + weight *= geometry.GetPeriodicEdgeWeight(iEdge); 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 10e0bf538cbd..40eb84479b0e 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -791,8 +791,9 @@ class CFVMFlowSolverBase : public CSolver { jPoint_UndLapl[iPoint] = msw ? 1 : 0; /*--- Loop over the neighbors of point i. ---*/ - for (auto jPoint : geometry->nodes->GetPoints(iPoint)) - { + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint); /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ @@ -808,7 +809,7 @@ class CFVMFlowSolverBase : public CSolver { } else { /*--- Jameson dissipation sensor, add variable difference and variable sum. ---*/ const auto weight = config->GetnMarker_Periodic() > 2 ? - geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) : 1.0; + geometry->GetPeriodicEdgeWeight(iEdge) : 1.0; iPoint_UndLapl[iPoint] += weight * (sensVar_j - sensVar_i); jPoint_UndLapl[iPoint] += weight * (sensVar_j + sensVar_i); } diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 85a290d6a4a5..304c665d0e25 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2397,7 +2397,9 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *co nodes->SetUnd_Lapl(iPoint, iVar, 0.0); /*--- Loop over the neighbors of point i. ---*/ - for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint); @@ -2405,7 +2407,7 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *co if (boundary_i && !boundary_j) continue; const auto weight = config->GetnMarker_Periodic() > 2 ? - geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) : 1.0; + geometry->GetPeriodicEdgeWeight(iEdge) : 1.0; /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index b2a3094283a0..4254d2318768 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -385,8 +385,10 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, 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); + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); + periodicNeighborCount(iPoint) += geometry->GetPeriodicEdgeWeight(iEdge); + } } END_SU2_OMP_FOR } @@ -587,13 +589,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, } nNeighbor = 0; - for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); /*--- Check if this neighbor lies on the periodic face so that we avoid double counting neighbors on both sides. If not, increment the count of neighbors for the donor. ---*/ - if (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) + if (geometry->GetPeriodicEdgeWeight(iEdge) == 1.0) nNeighbor++; } @@ -716,12 +719,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) Und_Lapl[iVar] = 0.0; - for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); /*--- Avoid periodic boundary points so that we do not duplicate edges on both sides of the periodic BC. ---*/ - if (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) { + if (geometry->GetPeriodicEdgeWeight(iEdge) == 1.0) { /*--- Solution differences ---*/ @@ -787,12 +792,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, on both sides of the periodic face. ---*/ Sensor_i = 0; Sensor_j = 0; - for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); /*--- Avoid halos and boundary points so that we don't duplicate edges on both sides of the periodic BC. ---*/ - if (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) < 1.0) continue; + if (geometry->GetPeriodicEdgeWeight(iEdge) < 1.0) continue; /*--- Use density instead of pressure for incomp. flows. ---*/ @@ -966,12 +973,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, r11 = 0.0; r12 = 0.0; r22 = 0.0; r13 = 0.0; r23_a = 0.0; r23_b = 0.0; r33 = 0.0; - for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); /*--- Avoid periodic boundary points so that we do not duplicate edges on both sides of the periodic BC. ---*/ - if (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) { + if (geometry->GetPeriodicEdgeWeight(iEdge) == 1.0) { /*--- Get coordinates for the neighbor point. ---*/ @@ -2460,7 +2469,9 @@ void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config) base_nodes->SetUnd_Lapl(iPoint, iVar, 0.0); /*--- Loop over the neighbors of point i. ---*/ - for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint); @@ -2471,7 +2482,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); + if (config->GetnMarker_Periodic() > 2) delta *= geometry->GetPeriodicEdgeWeight(iEdge); base_nodes->AddUnd_Lapl(iPoint, iVar, delta); } } diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index 671a8e1e262f..0e6bf2cebd13 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -271,6 +271,39 @@ TEST_CASE("Intersecting periodic stencils", "[Periodic][Gradients]") { } } +TEST_CASE("Periodic edge weight cache rebuild", "[Periodic][Geometry]") { + 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.InitConfig(); + field.InitGeometry(); + auto& geometry = *field.geometry; + geometry.SetPeriodicEdgeWeights(*field.config); + bool checked = false; + for (auto iEdge = 0ul; iEdge < geometry.GetnEdge(); ++iEdge) { + if (geometry.GetPeriodicEdgeWeight(iEdge) == 1.0) continue; + const auto iPoint = geometry.edges->GetNode(iEdge, 0); + for (auto iMarker = 0u; iMarker < geometry.GetnMarker(); ++iMarker) { + if (field.config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) continue; + const auto iVertex = geometry.nodes->GetVertex(iPoint, iMarker); + if (iVertex < 0) continue; + CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 0.5); + geometry.nodes->SetVertex(iPoint, -1, iMarker); + geometry.SetPeriodicEdgeWeights(*field.config); + CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 1.0); + geometry.nodes->SetVertex(iPoint, iVertex, iMarker); + geometry.SetPeriodicEdgeWeights(*field.config); + CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 0.5); + checked = true; + break; + } + break; + } + REQUIRE(checked); +} + TEST_CASE("Coarse periodic boundary flags", "[Periodic]") { UnitQuadTestCase field; const auto start = field.config_options.find("MARKER_HEATFLUX="); From 91393ce35df090e039d741e29d9d8da792230ce2 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 18:36:31 +0200 Subject: [PATCH 08/16] Share the principal periodic rotation angle calculation --- Common/include/toolboxes/geometry_toolbox.hpp | 14 ++++++++++++++ Common/src/geometry/CGeometry.cpp | 12 ++++++++---- UnitTests/SU2_CFD/gradients.cpp | 13 +++++++++++++ 3 files changed, 35 insertions(+), 4 deletions(-) diff --git a/Common/include/toolboxes/geometry_toolbox.hpp b/Common/include/toolboxes/geometry_toolbox.hpp index e621ea65b7b2..6e9e21da7398 100644 --- a/Common/include/toolboxes/geometry_toolbox.hpp +++ b/Common/include/toolboxes/geometry_toolbox.hpp @@ -206,6 +206,20 @@ inline void RotationMatrix(Scalar psi, Matrix& mat) { mat[1][1] = cosPsi; } +/*! \return Principal rotation angle in radians, using only the z angle in 2D. */ +template +inline Scalar RotationAngle(unsigned short nDim, const Scalar* angles) { + Scalar cosAngle; + if (nDim == 2) { + cosAngle = cos(angles[2]); + } else { + Scalar rotation[3][3]; + RotationMatrix(angles[0], angles[1], angles[2], rotation); + cosAngle = 0.5 * (rotation[0][0] + rotation[1][1] + rotation[2][2] - 1.0); + } + return acos(std::max(Scalar(-1.0), std::min(Scalar(1.0), cosAngle))); +} + /*! \brief Apply a rotation matrix (R) about origin (O) to a point at * distance (d) from it to obtain new coordinate (c). */ template diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 400cb6e99beb..b1c5676174a7 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -1172,10 +1172,14 @@ void CGeometry::SetPeriodicEdgeWeights(const CConfig& config) { * 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); + const auto angle = SU2_TYPE::GetValue(GeometryToolbox::RotationAngle(nDim, angles)); + constexpr passivedouble angleTol = 1e-4; + const auto nCopy = angle > angleTol ? static_cast(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 " + tag + ").", + CURRENT_FUNCTION); + weight *= 1.0 / nCopy; } else { weight *= 0.5; } diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index 0e6bf2cebd13..2e39d90a19d7 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -271,6 +271,19 @@ TEST_CASE("Intersecting periodic stencils", "[Periodic][Gradients]") { } } +TEST_CASE("Principal periodic rotation angle", "[Periodic][Geometry]") { + const su2double zero[3] = {}; + CHECK(GeometryToolbox::RotationAngle(3, zero) == 0.0); + const su2double quarterTurn[3] = {0, 0, -PI_NUMBER / 2}; + CHECK(GeometryToolbox::RotationAngle(2, quarterTurn) == Approx(PI_NUMBER / 2)); + CHECK(GeometryToolbox::RotationAngle(3, quarterTurn) == Approx(PI_NUMBER / 2)); + const su2double composed[3] = {PI_NUMBER / 2, PI_NUMBER / 2, 0}; + CHECK(GeometryToolbox::RotationAngle(2, composed) == 0.0); + CHECK(GeometryToolbox::RotationAngle(3, composed) == Approx(2 * PI_NUMBER / 3)); + const su2double halfTurn[3] = {PI_NUMBER, 0, 0}; + CHECK(GeometryToolbox::RotationAngle(3, halfTurn) == Approx(PI_NUMBER)); +} + TEST_CASE("Periodic edge weight cache rebuild", "[Periodic][Geometry]") { UnitQuadTestCase field; const auto start = field.config_options.find("MARKER_HEATFLUX="); From 71fc3863e1ebc2ddf366a1e18ea61d01dad30099 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 18:38:31 +0200 Subject: [PATCH 09/16] Detect periodic pair intersections within each mesh zone --- Common/include/geometry/CGeometry.hpp | 6 ++- Common/src/geometry/CGeometry.cpp | 35 ++++++++++++- .../computeGradientsLeastSquares.hpp | 2 +- .../include/solvers/CFVMFlowSolverBase.hpp | 2 +- SU2_CFD/src/solvers/CEulerSolver.cpp | 2 +- SU2_CFD/src/solvers/CSolver.cpp | 6 +-- .../check_periodic_intersection_rejection.py | 18 +++++++ UnitTests/SU2_CFD/gradients.cpp | 50 +++++++++++++++++-- UnitTests/meson.build | 3 ++ 9 files changed, 112 insertions(+), 12 deletions(-) create mode 100644 UnitTests/SU2_CFD/check_periodic_intersection_rejection.py diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 7ab090845020..bc909efa1cb4 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -319,6 +319,7 @@ class CGeometry { PassiveRequest* reqP_P2PSend{nullptr}; /*!< \brief Data structure for point-to-point send requests. */ PassiveRequest* reqP_P2PRecv{nullptr}; /*!< \brief Data structure for point-to-point recv requests. */ + bool periodicPairIntersections{false}; /*!< Whether distinct periodic pairs intersect in this zone. */ vector periodicEdgeWeight; /*!< Fraction of each edge stencil owned by this periodic copy. */ /*--- Data structures for periodic communications. ---*/ @@ -523,7 +524,10 @@ class CGeometry { /*! \brief Rebuild after changing connectivity, periodic membership or transforms; not needed for coordinate motion. */ - void SetPeriodicEdgeWeights(const CConfig& config); + void PreprocessPeriodicGeometry(const CConfig& config); + + /*! \brief Zone-wide decision, identical on all MPI ranks, for intersecting periodic stencils. */ + inline bool HasPeriodicPairIntersections() const { return periodicPairIntersections; } /*! \brief Fraction of an edge stencil owned by this periodic copy. */ inline passivedouble GetPeriodicEdgeWeight(unsigned long iEdge) const { diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index b1c5676174a7..8af884c37fda 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -1155,7 +1155,38 @@ void CGeometry::MatchPeriodic(const CConfig* config, unsigned short val_periodic delete[] Buffer_Recv_Marker; } -void CGeometry::SetPeriodicEdgeWeights(const CConfig& config) { +void CGeometry::PreprocessPeriodicGeometry(const CConfig& config) { + const auto nPairs = config.GetnMarker_Periodic() / 2; + if (nPairs == 0) { + periodicPairIntersections = false; + periodicEdgeWeight.clear(); + return; + } + + vector firstPair(nPoint, 0); + vector rotating(nPoint, false); + unsigned short local[2] = {}, global[2] = {}; + for (auto iMarker = 0u; iMarker < nMarker; ++iMarker) { + if (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) continue; + const auto pair = (config.GetMarker_All_PerBound(iMarker) - 1) % nPairs + 1; + const auto rotation = + GeometryToolbox::HasRotation(config.GetPeriodicRotAngles(config.GetMarker_All_TagBound(iMarker))); + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; ++iVertex) { + const auto iPoint = vertex[iMarker][iVertex]->GetNode(); + if (nodes->GetVertex(iPoint, iMarker) < 0) continue; + if (firstPair[iPoint] != 0 && firstPair[iPoint] != pair) { + local[0] = 1; + if (rotation || rotating[iPoint]) local[1] = 1; + } + if (firstPair[iPoint] == 0) firstPair[iPoint] = pair; + rotating[iPoint] = rotating[iPoint] || rotation; + } + } + SU2_MPI::Allreduce(local, global, 2, MPI_UNSIGNED_SHORT, MPI_MAX, SU2_MPI::GetComm()); + periodicPairIntersections = global[0] != 0; + if (global[1]) + SU2_MPI::Error("Intersecting periodic pairs involving rotation are not supported yet.", CURRENT_FUNCTION); + periodicEdgeWeight.assign(nEdge, 1.0); for (auto iEdge = 0ul; iEdge < nEdge; ++iEdge) { const auto iPoint = edges->GetNode(iEdge, 0), jPoint = edges->GetNode(iEdge, 1); @@ -1189,7 +1220,7 @@ void CGeometry::SetPeriodicEdgeWeights(const CConfig& config) { } void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { - geometry->SetPeriodicEdgeWeights(*config); + geometry->PreprocessPeriodicGeometry(*config); /*--- We start with the send and receive lists already available in the form of stored periodic point-donor pairs. We will loop through diff --git a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp index eb259db87132..8be32fa02974 100644 --- a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp +++ b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp @@ -254,7 +254,7 @@ void computeGradientsLeastSquares(CSolver* solver, if (weight > 0.0) { weight = 1.0 / weight; - if (config.GetnMarker_Periodic() > 2) + if (geometry.HasPeriodicPairIntersections()) weight *= geometry.GetPeriodicEdgeWeight(iEdge); for (size_t iDim = 0; iDim < nDim; ++iDim) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 40eb84479b0e..73ca86c3da83 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -808,7 +808,7 @@ 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. ---*/ - const auto weight = config->GetnMarker_Periodic() > 2 ? + const auto weight = geometry->HasPeriodicPairIntersections() ? geometry->GetPeriodicEdgeWeight(iEdge) : 1.0; iPoint_UndLapl[iPoint] += weight * (sensVar_j - sensVar_i); jPoint_UndLapl[iPoint] += weight * (sensVar_j + sensVar_i); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 304c665d0e25..b27825ebac8b 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2406,7 +2406,7 @@ 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 ? + const auto weight = geometry->HasPeriodicPairIntersections() ? geometry->GetPeriodicEdgeWeight(iEdge) : 1.0; /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 4254d2318768..01ce74d9fc4e 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -379,7 +379,7 @@ 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; + const bool intersectingPairs = geometry->HasPeriodicPairIntersections(); 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) @@ -1305,7 +1305,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, break; case PERIODIC_NEIGHBORS: - if (config->GetnMarker_Periodic() > 2) { + if (geometry->HasPeriodicPairIntersections()) { periodicNeighborCount(iPoint) += bufDRecv[buf_offset]; geometry->nodes->SetnNeighbor(iPoint, SU2_TYPE::Int(periodicNeighborCount(iPoint) + 0.5)); } else { @@ -2482,7 +2482,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(iEdge); + if (geometry->HasPeriodicPairIntersections()) delta *= geometry->GetPeriodicEdgeWeight(iEdge); base_nodes->AddUnd_Lapl(iPoint, iVar, delta); } } diff --git a/UnitTests/SU2_CFD/check_periodic_intersection_rejection.py b/UnitTests/SU2_CFD/check_periodic_intersection_rejection.py new file mode 100644 index 000000000000..e1ac8eb0f751 --- /dev/null +++ b/UnitTests/SU2_CFD/check_periodic_intersection_rejection.py @@ -0,0 +1,18 @@ +"""Run the MPI-aborting geometry check in a separate process.""" + +import subprocess +import sys + +result = subprocess.run( + sys.argv[1:] + ["Invalid rotational intersection"], + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + text=True, + timeout=60, +) +print(result.stdout) +assert result.returncode != 0, "Rotational intersection was accepted" +assert ( + "Intersecting periodic pairs involving rotation are not supported yet." + in result.stdout +) diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index 2e39d90a19d7..c36ab717e221 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -200,6 +200,7 @@ TEST_CASE("Intersecting periodic stencils", "[Periodic][Gradients]") { } for (auto pair = 1u; pair <= nPairs; ++pair) field.geometry->MatchPeriodic(field.config.get(), pair); field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + CHECK(field.geometry->HasPeriodicPairIntersections() == (field.config->GetnMarker_Periodic() > 2)); field.InitSolver(); auto* solver = field.solver[FLOW_SOL]; auto* nodes = solver->GetNodes(); @@ -271,6 +272,48 @@ TEST_CASE("Intersecting periodic stencils", "[Periodic][Gradients]") { } } +TEST_CASE("Zone-local periodic pair intersections", "[Periodic][Geometry]") { + const bool absentPair = 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, + absentPair ? "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n" : "MARKER_CUSTOM= (z_plus,z_minus)\n"); + field.AddOption(absentPair ? "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0, " + "other_minus,other_plus, 0,0,0, 0,0,0, 0,1,0)" + : "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)"); + field.InitConfig(); + field.InitGeometry(); + REQUIRE(field.config->GetnMarker_Periodic() == 4); + if (!absentPair && SU2_MPI::GetRank() != 0) { + /*--- Only rank zero has a local intersection; every rank must select the same communication path. ---*/ + for (auto marker = 0u; marker < field.geometry->GetnMarker(); ++marker) { + if (field.config->GetMarker_All_PerBound(marker) != 2 && field.config->GetMarker_All_PerBound(marker) != 4) + continue; + for (auto v = 0ul; v < field.geometry->GetnVertex(marker); ++v) + field.geometry->nodes->SetVertex(field.geometry->vertex[marker][v]->GetNode(), -1, marker); + } + } + field.geometry->PreprocessPeriodicGeometry(*field.config); + CHECK(field.geometry->HasPeriodicPairIntersections() == !absentPair); +} + +TEST_CASE("Invalid rotational intersection", "[.PeriodicInvalid]") { + 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= (z_plus,z_minus)\n"); + field.AddOption( + "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,90, 1,0,0, " + "y_minus,y_plus, 0,0,0, 0,0,0, 0,1,0)"); + field.InitConfig(); + field.InitGeometry(); + field.geometry->PreprocessPeriodicGeometry(*field.config); + FAIL("Rotational intersection was accepted"); +} + TEST_CASE("Principal periodic rotation angle", "[Periodic][Geometry]") { const su2double zero[3] = {}; CHECK(GeometryToolbox::RotationAngle(3, zero) == 0.0); @@ -293,7 +336,7 @@ TEST_CASE("Periodic edge weight cache rebuild", "[Periodic][Geometry]") { field.InitConfig(); field.InitGeometry(); auto& geometry = *field.geometry; - geometry.SetPeriodicEdgeWeights(*field.config); + geometry.PreprocessPeriodicGeometry(*field.config); bool checked = false; for (auto iEdge = 0ul; iEdge < geometry.GetnEdge(); ++iEdge) { if (geometry.GetPeriodicEdgeWeight(iEdge) == 1.0) continue; @@ -304,10 +347,10 @@ TEST_CASE("Periodic edge weight cache rebuild", "[Periodic][Geometry]") { if (iVertex < 0) continue; CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 0.5); geometry.nodes->SetVertex(iPoint, -1, iMarker); - geometry.SetPeriodicEdgeWeights(*field.config); + geometry.PreprocessPeriodicGeometry(*field.config); CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 1.0); geometry.nodes->SetVertex(iPoint, iVertex, iMarker); - geometry.SetPeriodicEdgeWeights(*field.config); + geometry.PreprocessPeriodicGeometry(*field.config); CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 0.5); checked = true; break; @@ -364,6 +407,7 @@ TEST_CASE("One-cell periodic stencil", "[Periodic][Gradients]") { field.InitGeometry(true); field.geometry->MatchPeriodic(field.config.get(), 1); field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + CHECK(field.geometry->HasPeriodicPairIntersections() == (field.config->GetnMarker_Periodic() > 2)); field.InitSolver(); auto* solver = field.solver[FLOW_SOL]; auto* nodes = solver->GetNodes(); diff --git a/UnitTests/meson.build b/UnitTests/meson.build index c6216606338f..8f67b4ef016f 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -53,6 +53,9 @@ if get_option('enable-tests') cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args] ) test('Catch2 test driver', test_driver) + test('Reject intersecting rotational periodic pairs', python, + args: files('SU2_CFD/check_periodic_intersection_rejection.py') + [test_driver.full_path()], + timeout: 90) test('NEMO convergence field migration', python, args: files('SU2_CFD/output/check_nemo_convergence_migration.py') + [test_driver.full_path()], timeout: 300) From c7dc21d2d5f184d2fde86adda5284ec1b1ed7752 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 18:40:40 +0200 Subject: [PATCH 10/16] Share least-squares normal matrix storage and packing --- .../computeGradientsLeastSquares.hpp | 35 +++++++++-- SU2_CFD/src/solvers/CSolver.cpp | 24 ++----- UnitTests/SU2_CFD/gradients.cpp | 62 +++++++++++++++++++ 3 files changed, 96 insertions(+), 25 deletions(-) diff --git a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp index 8be32fa02974..7df44db54b9a 100644 --- a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp +++ b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp @@ -32,6 +32,32 @@ namespace detail { +/*! \brief Visit the upper triangle and the extra xz entry used by the 3D LS factorization. */ +template +FORCEINLINE void forEachLeastSquaresNormalEntry(Function function) { + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = iDim; jDim < nDim; ++jDim) + function(iDim, jDim, iDim, jDim); + if (nDim == 3) function(2, 1, 0, 2); +} + +/*! \brief Rotate the symmetric normal matrix and pack it in the gradient storage convention. */ +template +FORCEINLINE void rotateLeastSquaresNormalMatrix(size_t iPoint, const MatrixType& matrix, + const su2double (&rotation)[nDim][nDim], su2double* packed) { + su2double rotated[nDim][nDim] = {}; + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = iDim; jDim < nDim; ++jDim) + for (size_t kDim = 0; kDim < nDim; ++kDim) + for (size_t lDim = 0; lDim < nDim; ++lDim) + rotated[iDim][jDim] += rotation[iDim][kDim] * matrix(iPoint, min(kDim,lDim), max(kDim,lDim)) * + rotation[jDim][lDim]; + std::fill_n(packed, nDim * nDim, su2double(0)); + forEachLeastSquaresNormalEntry([&](size_t iDim, size_t jDim, size_t row, size_t col) { + packed[iDim * nDim + jDim] = rotated[row][col]; + }); +} + /*! * \brief Prepare Smatrix for 2D. * \ingroup FvmAlgos @@ -257,12 +283,9 @@ void computeGradientsLeastSquares(CSolver* solver, if (geometry.HasPeriodicPairIntersections()) weight *= geometry.GetPeriodicEdgeWeight(iEdge); - for (size_t iDim = 0; iDim < nDim; ++iDim) - for (size_t jDim = iDim; jDim < nDim; ++jDim) - Rmatrix(iPoint,iDim,jDim) += dist_ij[iDim]*dist_ij[jDim]*weight; - - if (nDim == 3) - Rmatrix(iPoint,2,1) += dist_ij[0]*dist_ij[nDim-1]*weight; + forEachLeastSquaresNormalEntry([&](size_t iDim, size_t jDim, size_t row, size_t col) { + Rmatrix(iPoint,iDim,jDim) += dist_ij[row]*dist_ij[col]*weight; + }); /*--- Entries of c:= transpose(A)*b ---*/ diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 01ce74d9fc4e..4d8547a4abd3 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -898,25 +898,11 @@ 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 iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) - for (auto jDim = 0u; jDim < nDim; ++jDim) - 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[iCoordinate][kDim] : rotMatrix3D[iCoordinate][kDim]; - const auto qjl = nDim == 2 ? rotMatrix2D[jDim][lDim] : rotMatrix3D[jDim][lDim]; - rotated[iCoordinate][jDim] += qik * matrix[kDim][lDim] * qjl; - } - for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) - for (auto jDim = 0u; jDim < nDim; ++jDim) - bufDSend[buf_offset++] = iCoordinate <= jDim ? rotated[iCoordinate][jDim] : - (nDim == 3 && iCoordinate == 2 && jDim == 1 ? rotated[0][2] : su2double(0)); + if (nDim == 2) + detail::rotateLeastSquaresNormalMatrix(iPoint, base_nodes->GetRmatrix(), rotMatrix2D, &bufDSend[buf_offset]); + else + detail::rotateLeastSquaresNormalMatrix(iPoint, base_nodes->GetRmatrix(), rotMatrix3D, &bufDSend[buf_offset]); + buf_offset += nDim * nDim; for (auto iField = 0u; iField < ICOUNT; ++iField) Rotate(zeros, gradient[iPoint][iField], rotBlock[iField]); if (rotate_periodic) { diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index c36ab717e221..2fa552e20702 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -272,6 +272,68 @@ TEST_CASE("Intersecting periodic stencils", "[Periodic][Gradients]") { } } +template +void checkRotatedLeastSquaresNormalMatrix() { + C3DDoubleMatrix matrix(1, nDim, nDim, 0.0); + const su2double vectors[][3] = {{1, 2, 3}, {2, -1, 1}, {-1, 1, 2}}; + su2double rotation[nDim][nDim] = {}; + if constexpr (nDim == 2) + GeometryToolbox::RotationMatrix(su2double(PI_NUMBER / 2), rotation); + else + GeometryToolbox::RotationMatrix(su2double(0), su2double(0), su2double(PI_NUMBER / 2), rotation); + su2double expected[nDim][nDim] = {}; + for (const auto& vector : vectors) { + detail::forEachLeastSquaresNormalEntry( + [&](size_t i, size_t j, size_t row, size_t col) { matrix(0, i, j) += vector[row] * vector[col]; }); + su2double rotated[nDim] = {}, zeros[nDim] = {}; + GeometryToolbox::Rotate(rotation, zeros, vector, rotated); + for (auto i = 0u; i < nDim; ++i) + for (auto j = 0u; j < nDim; ++j) expected[i][j] += rotated[i] * rotated[j]; + } + su2double packed[nDim * nDim]; + detail::rotateLeastSquaresNormalMatrix(0, matrix, rotation, packed); + for (auto i = 0u; i < nDim; ++i) + for (auto j = 0u; j < nDim; ++j) { + const auto ref = i <= j ? expected[i][j] : (nDim == 3 && i == 2 && j == 1 ? expected[0][nDim - 1] : 0.0); + CHECK(packed[i * nDim + j] == Approx(ref).margin(1e-12)); + } +} + +TEST_CASE("Rotated least-squares normal matrix", "[Periodic][Gradients]") { + checkRotatedLeastSquaresNormalMatrix<2>(); + checkRotatedLeastSquaresNormalMatrix<3>(); +} + +TEST_CASE("2D periodic linear gradient", "[Periodic][Gradients]") { + const bool weighted = GENERATE(false, true); + UnitQuadTestCase field; + field.SetOption("MESH_FORMAT= RECTANGLE"); + field.SetOption("MARKER_CUSTOM= (z_plus,z_minus)"); + field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); + 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) + for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) + nodes->SetPrimitive(i, v, 3 + field.geometry->nodes->GetCoord(i, 1)); + solver->SetRotatePeriodic(false); + SU2_OMP_PARALLEL { + computeGradientsLeastSquares(solver, MPI_QUANTITIES::PRIMITIVE_GRADIENT, + weighted ? PERIODIC_PRIM_LS : PERIODIC_PRIM_ULS, *field.geometry, *field.config, + weighted, nodes->GetPrimitive(), 0, solver->GetnPrimVarGrad(), -1, + nodes->GetGradient_Primitive(), nodes->GetRmatrix()); + } + END_SU2_OMP_PARALLEL + for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) { + CHECK(nodes->GetGradient_Primitive(i, 0, 0) == Approx(0).margin(1e-12)); + CHECK(nodes->GetGradient_Primitive(i, 0, 1) == Approx(1).margin(1e-12)); + } +} + TEST_CASE("Zone-local periodic pair intersections", "[Periodic][Geometry]") { const bool absentPair = GENERATE(false, true); UnitQuadTestCase field; From 23f4300bb5d259197a92e11ccf1053d38cbf0271 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 18:40:41 +0200 Subject: [PATCH 11/16] Remove obsolete periodic unsigned-short buffers --- Common/include/geometry/CGeometry.hpp | 2 -- Common/src/geometry/CGeometry.cpp | 23 ----------------------- 2 files changed, 25 deletions(-) diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index bc909efa1cb4..a12ec045164e 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -349,8 +349,6 @@ class CGeometry { 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. */ - 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. */ SU2_MPI::Request* req_PeriodicRecv{nullptr}; /*!< \brief Data structure for periodic recv requests. */ diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 8af884c37fda..446b98691677 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -127,9 +127,6 @@ CGeometry::~CGeometry() { delete[] bufD_PeriodicRecv; delete[] bufD_PeriodicSend; - delete[] bufS_PeriodicRecv; - delete[] bufS_PeriodicSend; - delete[] req_PeriodicSend; delete[] req_PeriodicRecv; @@ -1379,9 +1376,6 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { bufD_PeriodicSend = nullptr; bufD_PeriodicRecv = nullptr; - bufS_PeriodicSend = nullptr; - bufS_PeriodicRecv = nullptr; - /*--- Allocate memory for the MPI requests if we need to communicate. ---*/ if (nPeriodicSend > 0) { @@ -1584,12 +1578,6 @@ void CGeometry::AllocatePeriodicComms(unsigned short countPerPeriodicPoint) { delete[] bufD_PeriodicRecv; bufD_PeriodicRecv = new su2double[nRecv](); - - delete[] bufS_PeriodicSend; - bufS_PeriodicSend = new unsigned short[nSend](); - - delete[] bufS_PeriodicRecv; - bufS_PeriodicRecv = new unsigned short[nRecv](); } END_SU2_OMP_SAFE_GLOBAL_ACCESS } @@ -1631,10 +1619,6 @@ void CGeometry::PostPeriodicRecvs(CGeometry* geometry, const CConfig* config, CO SU2_MPI::Irecv(&(static_cast(bufD_PeriodicRecv)[offset]), count, MPI_DOUBLE, source, tag, SU2_MPI::GetComm(), &(req_PeriodicRecv[iRecv])); break; - case COMM_TYPE::UNSIGNED_SHORT: - SU2_MPI::Irecv(&(static_cast(bufS_PeriodicRecv)[offset]), count, MPI_UNSIGNED_SHORT, source, - tag, SU2_MPI::GetComm(), &(req_PeriodicRecv[iRecv])); - break; default: SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION); break; @@ -1678,10 +1662,6 @@ void CGeometry::PostPeriodicSends(CGeometry* geometry, const CConfig* config, CO SU2_MPI::Isend(&(static_cast(bufD_PeriodicSend)[offset]), count, MPI_DOUBLE, dest, tag, SU2_MPI::GetComm(), &(req_PeriodicSend[val_iSend])); break; - case COMM_TYPE::UNSIGNED_SHORT: - SU2_MPI::Isend(&(static_cast(bufS_PeriodicSend)[offset]), count, MPI_UNSIGNED_SHORT, dest, tag, - SU2_MPI::GetComm(), &(req_PeriodicSend[val_iSend])); - break; default: SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION); break; @@ -1702,9 +1682,6 @@ void CGeometry::PostPeriodicSends(CGeometry* geometry, const CConfig* config, CO case COMM_TYPE::DOUBLE: parallelCopy(myFinal - myStart, &bufD_PeriodicSend[myStart], &bufD_PeriodicRecv[iRecv]); break; - case COMM_TYPE::UNSIGNED_SHORT: - parallelCopy(myFinal - myStart, &bufS_PeriodicSend[myStart], &bufS_PeriodicRecv[iRecv]); - break; default: SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION); break; From deb3cb87396e5f5dae09acf9e1b4e3b3758f8cfe Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 18:40:41 +0200 Subject: [PATCH 12/16] Rotate periodic Laplacian directly into the send buffer --- SU2_CFD/src/solvers/CSolver.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 4d8547a4abd3..245137b0bfdb 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -703,12 +703,10 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, case PERIODIC_LAPLACIAN: if (intersectingPairs) { + const auto* undLapl = base_nodes->GetUndivided_Laplacian(iPoint); 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 iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) - bufDSend[buf_offset + 1 + iCoordinate] = Und_Lapl[1 + iCoordinate]; + bufDSend[buf_offset + iField] = undLapl[iField]; + if (rotate_periodic) Rotate(zeros, &undLapl[1], &bufDSend[buf_offset + 1]); break; } From cc4c4cd38eacbf7ffe7ae66a807cce3aecf9dedc Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 18:54:03 +0200 Subject: [PATCH 13/16] Keep cache rebuild regression collectives aligned across MPI ranks --- UnitTests/SU2_CFD/gradients.cpp | 35 ++++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index 2fa552e20702..95b37db5cf9e 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -396,10 +396,12 @@ TEST_CASE("Periodic edge weight cache rebuild", "[Periodic][Geometry]") { 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.InitConfig(); - field.InitGeometry(); + field.InitGeometry(true); auto& geometry = *field.geometry; geometry.PreprocessPeriodicGeometry(*field.config); - bool checked = false; + unsigned long edge = geometry.GetnEdge(), point = 0; + unsigned short marker = 0; + long vertex = -1; for (auto iEdge = 0ul; iEdge < geometry.GetnEdge(); ++iEdge) { if (geometry.GetPeriodicEdgeWeight(iEdge) == 1.0) continue; const auto iPoint = geometry.edges->GetNode(iEdge, 0); @@ -407,19 +409,30 @@ TEST_CASE("Periodic edge weight cache rebuild", "[Periodic][Geometry]") { if (field.config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) continue; const auto iVertex = geometry.nodes->GetVertex(iPoint, iMarker); if (iVertex < 0) continue; - CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 0.5); - geometry.nodes->SetVertex(iPoint, -1, iMarker); - geometry.PreprocessPeriodicGeometry(*field.config); - CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 1.0); - geometry.nodes->SetVertex(iPoint, iVertex, iMarker); - geometry.PreprocessPeriodicGeometry(*field.config); - CHECK(geometry.GetPeriodicEdgeWeight(iEdge) == 0.5); - checked = true; + edge = iEdge; + point = iPoint; + marker = iMarker; + vertex = iVertex; break; } break; } - REQUIRE(checked); + const bool checked = edge < geometry.GetnEdge(); + if (checked) { + CHECK(geometry.GetPeriodicEdgeWeight(edge) == 0.5); + geometry.nodes->SetVertex(point, -1, marker); + } + /*--- All ranks rebuild together, including those without a local candidate edge. ---*/ + geometry.PreprocessPeriodicGeometry(*field.config); + if (checked) { + CHECK(geometry.GetPeriodicEdgeWeight(edge) == 1.0); + geometry.nodes->SetVertex(point, vertex, marker); + } + geometry.PreprocessPeriodicGeometry(*field.config); + if (checked) CHECK(geometry.GetPeriodicEdgeWeight(edge) == 0.5); + unsigned long local = checked, global = 0; + SU2_MPI::Allreduce(&local, &global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(global > 0); } TEST_CASE("Coarse periodic boundary flags", "[Periodic]") { From fbb7d3a86276a1bb2602fcaa351d9ebc1fc45e31 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 20:38:14 +0200 Subject: [PATCH 14/16] Promote periodic least-squares buffer stride before multiplying --- SU2_CFD/src/solvers/CSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 245137b0bfdb..90353cbe9745 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -900,7 +900,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, detail::rotateLeastSquaresNormalMatrix(iPoint, base_nodes->GetRmatrix(), rotMatrix2D, &bufDSend[buf_offset]); else detail::rotateLeastSquaresNormalMatrix(iPoint, base_nodes->GetRmatrix(), rotMatrix3D, &bufDSend[buf_offset]); - buf_offset += nDim * nDim; + buf_offset += static_cast(nDim) * nDim; for (auto iField = 0u; iField < ICOUNT; ++iField) Rotate(zeros, gradient[iPoint][iField], rotBlock[iField]); if (rotate_periodic) { From 6a7fb5f94888e999f1ae905a4f5f8ad41c237dcb Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sun, 18 Oct 2026 20:38:59 +0200 Subject: [PATCH 15/16] Refresh serial and MPI Aachen references after periodic stencil fixes --- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 4c8028e266dd..1e6b7b6395c8 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1331,7 +1331,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.688483, -8.466346, -6.035067, -6.425967, -5.820125, -4.597817, -5.528297, -5.301466, -3.819862, -5.242356, -5.738805, -3.630202, -2.213240, -2.857919, -0.576864] + Aachen_3D_restart.test_vals = [-7.688401, -8.466746, -6.034652, -6.427491, -5.818347, -4.597687, -5.530255, -5.303027, -3.819060, -5.242491, -5.755143, -3.629921, -2.213103, -2.859285, -0.576337] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 1edd0378ab2a..c211c1852ffa 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -957,7 +957,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.689117, -8.476204, -6.036107, -6.427402, -5.831576, -4.599024, -5.530512, -5.305687, -3.821789, -5.243255, -5.754831, -3.631603, -2.213934, -2.895728, -0.578545] + Aachen_3D_restart.test_vals = [-7.689029, -8.475341, -6.035704, -6.426165, -5.830215, -4.598893, -5.529338, -5.301964, -3.820993, -5.243202, -5.732910, -3.631319, -2.213792, -2.876325, -0.578057] Aachen_3D_restart.enabled_with_asan = False test_list.append(Aachen_3D_restart) From 977dc88b74216fb15ed60030855fdb0962f2f4b9 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 19 Oct 2026 14:57:29 +0200 Subject: [PATCH 16/16] Update regression references after the periodic stencil and coarse flag fixes Periodic cases change because of the stencil fixes. poiseuille_profile changes because coarse multigrid points at the inlet now inherit the viscous wall flag. The other cases change only by rounding of the rewritten loops. Co-Authored-By: Claude Opus 5.5 --- TestCases/hybrid_regression.py | 8 ++++---- TestCases/parallel_regression.py | 10 +++++----- TestCases/parallel_regression_AD.py | 2 +- TestCases/serial_regression.py | 8 ++++---- TestCases/serial_regression_AD.py | 2 +- TestCases/tutorials.py | 2 +- 6 files changed, 16 insertions(+), 16 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 21c63694f567..d298966939fb 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -149,7 +149,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.004080, -7.638116, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.004094, -7.640551, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.004276, -7.636719, -0.000000, 2.089953] test_list.append(poiseuille_profile) @@ -178,7 +178,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-1.461774, 5.883885, 0.630267, 0.020404, 100.000000] + rae2822_sst.test_vals = [-1.519695, 5.888685, 0.614691, 0.020767, 100.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -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.959252, -12.213981, -19.314222, -13.600248, -19.076384, -13.444335, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.961836, -12.211509, -19.315122, -13.598954, -19.081506, -13.442877, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -711,7 +711,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.455981, -1.541069] + bars_SST_2D.test_vals = [13.000000, -0.455976, -1.539981] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 1e6b7b6395c8..a0a2eeda708b 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -425,7 +425,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.003186, -7.623957, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.003181, -7.626400, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.007498, -7.226926, -0.000000, 2.089953] poiseuille_profile.tol = [0.001, 0.001, 1e-5, 1e-5, 1e-5] test_list.append(poiseuille_profile) @@ -965,7 +965,7 @@ def main(): turbmod_sa_neg_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_neg_rae2822.cfg_file = "turb_SA_NEG_RAE2822.cfg" turbmod_sa_neg_rae2822.test_iter = 10 - turbmod_sa_neg_rae2822.test_vals = [1.275994, 1.062588, -2.067822, 1.300662, 0.480881, 0.000000] + turbmod_sa_neg_rae2822.test_vals = [1.275906, 1.062762, -2.067335, 1.301619, 0.481422, 0.000000] turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.345593, 1.448310, 1.208721, -0.846597, 1.248410, 0.489117, 0.000000] test_list.append(turbmod_sa_neg_rae2822) @@ -1339,7 +1339,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.905382, -12.212100, -19.212143, -13.508268, -19.078783, -13.446389, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.903628, -12.211875, -19.212536, -13.508088, -19.079442, -13.440735, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -1447,7 +1447,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.457912, -1.541047] + bars_SST_2D.test_vals = [13.000000, -0.457907, -1.539958] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -1773,7 +1773,7 @@ def main(): mms_fvm_vortex.cfg_dir = "mms/fvm_euler" mms_fvm_vortex.cfg_file = "inv_mms_vortex.cfg" mms_fvm_vortex.test_iter = 10 - mms_fvm_vortex.test_vals = [-5.704300, -4.848072, 0.000000, 0.000000] + mms_fvm_vortex.test_vals = [-5.704011, -4.847997, 0.000000, 0.000000] mms_fvm_vortex.unsteady = True test_list.append(mms_fvm_vortex) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 250f96e66c12..ef8999cadd0e 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -315,7 +315,7 @@ def main(): discadj_fsi2.cfg_dir = "disc_adj_fsi/Airfoil_2d" discadj_fsi2.cfg_file = "config.cfg" discadj_fsi2.test_iter = 8 - discadj_fsi2.test_vals = [-3.824633, 1.979516, -3.863368, 0.295450, 3.839800] + discadj_fsi2.test_vals = [-3.824660, 1.979488, -3.863368, 0.295450, 3.839800] discadj_fsi2.test_vals_aarch64 = [-3.824870, 1.979160, -3.863368, 0.295450, 3.839800] discadj_fsi2.tol = 0.00001 test_list.append(discadj_fsi2) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index c211c1852ffa..e65b6586eae4 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -216,7 +216,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.003115, -7.626023, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.003109, -7.627917, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262299, -0.000000, 2.089953] #last 4 columns test_list.append(poiseuille_profile) @@ -966,7 +966,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.905642, -12.212181, -19.212554, -13.508053, -19.069991, -13.439494, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.935390, -12.212295, -19.248052, -13.539427, -19.086000, -13.455882, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -1086,7 +1086,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.456143, -1.541051] + bars_SST_2D.test_vals = [13.000000, -0.456139, -1.539962] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -1238,7 +1238,7 @@ def main(): mms_fvm_vortex.cfg_dir = "mms/fvm_euler" mms_fvm_vortex.cfg_file = "inv_mms_vortex.cfg" mms_fvm_vortex.test_iter = 10 - mms_fvm_vortex.test_vals = [-5.704300, -4.848072, 0.000000, 0.000000] + mms_fvm_vortex.test_vals = [-5.704011, -4.847997, 0.000000, 0.000000] mms_fvm_vortex.unsteady = True test_list.append(mms_fvm_vortex) diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 288f8a47efd3..44bfc4d0073d 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -200,7 +200,7 @@ def main(): discadj_trans_stator.cfg_dir = "disc_adj_turbomachinery/transonic_stator_2D" discadj_trans_stator.cfg_file = "transonic_stator.cfg" discadj_trans_stator.test_iter = 79 - discadj_trans_stator.test_vals = [79.000000, -7.308168, -9.891121, -10.038669, -13.368502] + discadj_trans_stator.test_vals = [79.000000, -7.308179, -9.891140, -10.038680, -13.368512] discadj_trans_stator.test_vals_aarch64 = [79.000000, -7.308167, -9.891117, -10.038669, -13.368501] test_list.append(discadj_trans_stator) diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index 22e9e55fdb03..f08b73f21d26 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -314,7 +314,7 @@ def main(): tutorial_unst_naca0012.cfg_dir = "../Tutorials/compressible_flow/Unsteady_NACA0012" tutorial_unst_naca0012.cfg_file = "unsteady_naca0012.cfg" tutorial_unst_naca0012.test_iter = 520 - tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.292186, 0.000000, 0.304255, 0.792893, 0.001907, 0.010888] + tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.296539, 0.000000, 0.308341, 0.797028, 0.003155, 0.015433] tutorial_unst_naca0012.test_vals_aarch64 = [520, 0, -5.292359, 0, 0.284720, 0.766329, 0.000954, 0.007565] tutorial_unst_naca0012.unsteady = True test_list.append(tutorial_unst_naca0012)