diff --git a/Common/include/toolboxes/geometry_toolbox.hpp b/Common/include/toolboxes/geometry_toolbox.hpp index 1087f940a36f..e621ea65b7b2 100644 --- a/Common/include/toolboxes/geometry_toolbox.hpp +++ b/Common/include/toolboxes/geometry_toolbox.hpp @@ -27,6 +27,7 @@ #pragma once #include +#include namespace GeometryToolbox { /// \addtogroup GeometryToolbox @@ -215,6 +216,31 @@ inline void Rotate(const Scalar R[][nDim], const Scalar* O, const Scalar* d, Sca } } +/*! \return Whether any of the three supplied rotation angles is nonzero. */ +template +inline bool HasRotation(const Scalar* angles) { + // Configured zero stays exact after degree conversion and negation; retain every nonzero rotation. + return angles[0] != 0.0 || angles[1] != 0.0 || angles[2] != 0.0; +} + +/*! \brief Rotate component bounds in place, enclosing the rotated box. */ +template +inline void RotateBox(const Scalar R[][nDim], Scalar* vMin, Scalar* vMax) { + using std::max; + using std::min; + Scalar rotMin[nDim] = {0.0}, rotMax[nDim] = {0.0}; + for (int iDim = 0; iDim < nDim; ++iDim) { + for (int jDim = 0; jDim < nDim; ++jDim) { + const Scalar fromMin = R[iDim][jDim] * vMin[jDim]; + const Scalar fromMax = R[iDim][jDim] * vMax[jDim]; + rotMin[iDim] += min(fromMin, fromMax); + rotMax[iDim] += max(fromMin, fromMax); + } + } + std::copy_n(rotMin, nDim, vMin); + std::copy_n(rotMax, nDim, vMax); +} + /*! \brief Tangent projection */ template inline void TangentProjection(Int nDim, const Mat& tensor, const Scalar* vector, Scalar* proj) { diff --git a/SU2_CFD/include/limiters/CLimiterDetails.hpp b/SU2_CFD/include/limiters/CLimiterDetails.hpp index 2fe298261006..554509f066fc 100644 --- a/SU2_CFD/include/limiters/CLimiterDetails.hpp +++ b/SU2_CFD/include/limiters/CLimiterDetails.hpp @@ -72,6 +72,16 @@ struct LimiterHelpers { FORCEINLINE static Type epsilon() {return std::numeric_limits::epsilon();} + /*! \brief MUSCL reconstruction increment using the displacement to the middle of an edge. */ + template + FORCEINLINE static Type reconstructionIncrement(Int nDim, const Type* halfEdge, const Type* gradient, + const Type& value_i, const Type& value_j, const Type& kappa) { + Type proj = 0.0; + for (Int iDim = 0; iDim < nDim; ++iDim) proj += halfEdge[iDim] * gradient[iDim]; + const Type cent = 0.5 * (value_j - value_i); + return umusclProjection(proj, cent, kappa); + } + FORCEINLINE static Type umusclProjection(const Type& grad_proj, const Type& delta, const Type& kappa) { /*-------------------------------------------------------------------*/ diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp index 7d8beafb03e0..3b19f241f3bc 100644 --- a/SU2_CFD/include/limiters/computeLimiters_impl.hpp +++ b/SU2_CFD/include/limiters/computeLimiters_impl.hpp @@ -109,11 +109,29 @@ void computeLimiters_impl(CSolver* solver, limiterDetails.preprocess(geometry, config, varBegin, varEnd, field); + /*--- With rotational periodicity, the first periodic comm. also brings the min/max + * projections over the edges of the periodic matches (stored after each other), + * because the limiters of the velocity cannot be compared across a rotation. ---*/ + + su2activematrix* periodicProj = nullptr; + + if (periodic && (kindPeriodicComm1 == PERIODIC_LIM_PRIM_1)) + periodicProj = solver->GetPeriodicProjections(geometry, config); + /*--- Initialize all min/max field values if we have * periodic comms. otherwise do it inside main loop. ---*/ if (periodic) { + if (periodicProj != nullptr) + { + SU2_OMP_FOR_STAT(chunkSize) + for (auto iPoint = 0ul; iPoint < periodicProj->rows(); ++iPoint) + for (auto iVar = 0ul; iVar < periodicProj->cols(); ++iVar) + (*periodicProj)(iPoint,iVar) = 0.0; + END_SU2_OMP_FOR + } + SU2_OMP_FOR_STAT(chunkSize) for (size_t iPoint = 0; iPoint < nPoint; ++iPoint) for (size_t iVar = varBegin; iVar < varEnd; ++iVar) @@ -165,6 +183,23 @@ void computeLimiters_impl(CSolver* solver, for (size_t iVar = varBegin; iVar < varEnd; ++iVar) projMax[iVar] = projMin[iVar] = 0.0; + if (periodicProj != nullptr && nodes->GetPeriodicBoundary(iPoint)) + { + /*--- Start from the min/max over the edges of the periodic matches. ---*/ + + const auto* projections = solver->GetPeriodicProjection(iPoint); + if (projections != nullptr) { + for (auto iVar = varBegin; iVar < varEnd; ++iVar) { + const auto& periodicMin = projections[iVar]; + const auto& periodicMax = projections[periodicProj->cols()/2 + iVar]; + AD::SetPreaccIn(periodicMin); + AD::SetPreaccIn(periodicMax); + projMin[iVar] = periodicMin; + projMax[iVar] = periodicMax; + } + } + } + /*--- Compute max/min projection and values over direct neighbors. ---*/ for (auto jPoint : geometry.nodes->GetPoints(iPoint)) { @@ -175,22 +210,16 @@ void computeLimiters_impl(CSolver* solver, /*--- Distance vector from iPoint to face (middle of the edge). ---*/ su2double dist_ij[nDim] = {0.0}; - - for(size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t iDim = 0; iDim < nDim; ++iDim) dist_ij[iDim] = 0.5 * (coord_j[iDim] - coord_i[iDim]); /*--- Project each variable, update min/max. ---*/ for(size_t iVar = varBegin; iVar < varEnd; ++iVar) { - su2double proj = 0.0; - - for(size_t iDim = 0; iDim < nDim; ++iDim) - proj += dist_ij[iDim] * gradient(iPoint,iVar,iDim); - AD::SetPreaccIn(field(jPoint,iVar)); - const su2double cent = 0.5 * (field(jPoint,iVar) - field(iPoint,iVar)); - proj = LimiterHelpers<>::umusclProjection(proj, cent, umusclKappa); + const su2double proj = LimiterHelpers<>::reconstructionIncrement(nDim, dist_ij, + gradient[iPoint][iVar], field(iPoint,iVar), field(jPoint,iVar), umusclKappa); projMax[iVar] = max(projMax[iVar], proj); projMin[iVar] = min(projMin[iVar], proj); @@ -224,7 +253,8 @@ void computeLimiters_impl(CSolver* solver, } END_SU2_OMP_FOR - /*--- Account for periodic effects, take the minimum limiter on each periodic pair. ---*/ + /*--- Account for periodic effects, take the minimum limiter on each periodic pair + * (except for the velocity with rotational periodicity). ---*/ if (periodic) { for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index cb44757c6280..3a46fc8c4db6 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -261,7 +261,7 @@ void CFVMFlowSolverBase::CommunicateInitialState(CGeometry* geometry, cons CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); } SetImplicitPeriodic(euler_implicit); - if (MGLevel == MESH_0) SetRotatePeriodic(true); + SetRotatePeriodic(true); /*--- Perform the MPI communication of the solution ---*/ diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 7a7b6c726021..76d883f2690b 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -37,6 +37,7 @@ #include #include #include +#include #include #include @@ -146,6 +147,8 @@ class CSolver { bool rotate_periodic; /*!< \brief Flag that controls whether the periodic solution needs to be rotated for the solver. */ bool implicit_periodic; /*!< \brief Flag that controls whether the implicit system should be treated by the periodic BC comms. */ + su2activematrix PeriodicProj; /*!< \brief Min/max reconstruction increments at periodic receive points (for limiters). */ + std::unordered_map PeriodicProjIndex; /*!< \brief Local point to compact projection row. */ bool dynamic_grid; /*!< \brief Flag that determines whether the grid is dynamic (moving or deforming + grid velocities). */ @@ -4231,6 +4234,22 @@ class CSolver { */ inline void SetRotatePeriodic(bool val_rotate_periodic) { rotate_periodic = val_rotate_periodic; } + /*! + * \brief Storage for the limiters with rotational periodicity: the min and max, over the edges of the periodic + * matches of each point, of the reconstruction increments (communicated with PERIODIC_LIM_PRIM_1). + * \param[in] geometry - Periodic receive points of this mesh level. + * \param[in] config - Definition of the particular problem. + * \return The matrix (unique periodic receive points x 2*nPrimVarGrad, min then max), nullptr without rotation. + */ + su2activematrix* GetPeriodicProjections(const CGeometry& geometry, const CConfig& config); + + /*! \brief Reconstruction increment bounds for a periodic receive point, nullptr for other points. */ + inline su2double* GetPeriodicProjection(unsigned long iPoint) { + const auto& indices = PeriodicProjIndex; + const auto it = indices.find(iPoint); + return it == indices.end() ? nullptr : PeriodicProj[it->second]; + } + /*! * \brief Retrieve the solver name for output purposes. * \returns Name of the solver. diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 56cee4b6ad52..fd02e9956918 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -262,7 +262,8 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, JCOUNT = nDim; break; case PERIODIC_LIM_PRIM_1: - COUNT_PER_POINT = nPrimVarGrad*2; + /*--- Min and max of the solution, min and max of the reconstruction increments. ---*/ + COUNT_PER_POINT = nPrimVarGrad*4; MPI_TYPE = COMM_TYPE::DOUBLE; ICOUNT = nPrimVarGrad; break; @@ -338,6 +339,30 @@ namespace PeriodicCommHelpers { } } +su2activematrix* CSolver::GetPeriodicProjections(const CGeometry& geometry, const CConfig& config) { + + if (!rotate_periodic) return nullptr; + + bool rotation = false; + for (auto iMarker = 0u; iMarker < config.GetnMarker_All(); iMarker++) { + if (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY) continue; + rotation |= GeometryToolbox::HasRotation(config.GetPeriodicRotAngles(config.GetMarker_All_TagBound(iMarker))); + } + if (!rotation) return nullptr; + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { + if (PeriodicProj.empty() && geometry.nPeriodicRecv > 0) { + for (auto iRecv = 0; iRecv < geometry.nPoint_PeriodicRecv[geometry.nPeriodicRecv]; ++iRecv) + PeriodicProjIndex.emplace(geometry.Local_Point_PeriodicRecv[iRecv], PeriodicProjIndex.size()); + PeriodicProj.resize(PeriodicProjIndex.size(), 2*nPrimVarGrad); + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + return &PeriodicProj; +} + void CSolver::InitiatePeriodicComms(CGeometry *geometry, const CConfig *config, unsigned short val_periodic_index, @@ -391,6 +416,13 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, else GeometryToolbox::Rotate(rotMatrix3D, origin, direction, rotated); }; + /*--- Rotates, in place, the range [vMin, vMax] of each component of a vector: each term + of the product with the rotation matrix takes its smallest and its largest value. ---*/ + auto RotateBox = [&](su2double* vMin, su2double* vMax) { + if (nDim == 2) GeometryToolbox::RotateBox(rotMatrix2D, vMin, vMax); + else GeometryToolbox::RotateBox(rotMatrix3D, vMin, vMax); + }; + string Marker_Tag; /*--- Set the size of the data packet and type depending on quantity. ---*/ @@ -421,6 +453,37 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, auto& limiter = PeriodicCommHelpers::selectLimiter(base_nodes, commType); auto& field = PeriodicCommHelpers::selectField(base_nodes, commType); + /*--- Updates the range [Sol_Min, Sol_Max] of each variable with the values of a neighbor, + the velocity of the neighbor (variables 1 to nDim) is rotated first if requested. ---*/ + auto UpdateMinMax = [&](const su2double* values, bool rotateVelocity) { + if (rotateVelocity) { + for (auto iField = 0u; iField < ICOUNT; iField++) rotPrim_j[iField] = values[iField]; + Rotate(zeros, &values[1], &rotPrim_j[1]); + values = rotPrim_j; + } + for (auto iField = 0u; iField < ICOUNT; iField++) { + Sol_Min[iField] = min(Sol_Min[iField], values[iField]); + Sol_Max[iField] = max(Sol_Max[iField], values[iField]); + } + }; + + /*--- Reconstruction increments of all variables from a point to the middle + of the edge to one of its neighbors, computed as in computeLimiters_impl. ---*/ + auto ReconstructionIncrements = [&](unsigned long point_i, unsigned long point_j, su2double* increments) { + const su2double kappa = config->GetMUSCL_Kappa_Flow(); + const auto* coord_i = geometry->nodes->GetCoord(point_i); + const auto* coord_j = geometry->nodes->GetCoord(point_j); + + su2double dist_ij[3] = {0.0}; + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + dist_ij[iCoordinate] = 0.5 * (coord_j[iCoordinate] - coord_i[iCoordinate]); + + for (auto iField = 0u; iField < ICOUNT; iField++) { + increments[iField] = LimiterHelpers<>::reconstructionIncrement(nDim, dist_ij, + gradient[point_i][iField], field(point_i, iField), field(point_j, iField), kappa); + } + }; + /*--- Load the specified quantity from the solver into the generic communication buffer in the geometry class. ---*/ @@ -479,6 +542,10 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, GeometryToolbox::RotationMatrix(Theta, Phi, Psi, rotMatrix3D); } + /*--- Whether the vector components of the solution are rotated for this marker. ---*/ + + const bool rotation = rotate_periodic && GeometryToolbox::HasRotation(angles); + /*--- Compute the offset in the recv buffer for this point. ---*/ buf_offset = (msg_offset + iSend)*COUNT_PER_POINT; @@ -553,7 +620,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, } } - /*--- Rotate the momentum columns of the Jacobian. ---*/ + /*--- Rotate the momentum rows and columns of the Jacobian, the residual of the + periodic match is Q*R(Q^T*U), so its Jacobian is Q*J*Q^T. First the rows. ---*/ if (rotate_periodic) { for (iVar = 0; iVar < nVar; iVar++) { @@ -569,6 +637,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, rotMatrix3D[2][2]*block(3, iVar); } } + + /*--- Then the columns, i.e. the momentum part of each row. ---*/ + + for (auto iRow = 0u; iRow < nVar; iRow++) { + su2double rotated[3] = {0.0}; + Rotate(zeros, &jacBlock[iRow][1], rotated); + for (auto jDim = 0u; jDim < nDim; jDim++) jacBlock[iRow][1+jDim] = rotated[jDim]; + } } /*--- Load the Jacobian terms into the buffer for sending. ---*/ @@ -936,11 +1012,15 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, Sol_Max[iVar] = base_nodes->GetSolution_Max()(iPoint, iVar); } + /*--- The velocity of each neighbour is rotated before taking the min and max of its + components (the rotated min and max vectors do not bound the rotated components). + The min/max stored so far (the value of the point, or what was accumulated from + other periodic pairs) are a range for each component, which is rotated as a box. ---*/ + + if (rotation) RotateBox(&Sol_Min[1], &Sol_Max[1]); + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { - for (iVar = 0; iVar < ICOUNT; iVar++) { - Sol_Min[iVar] = min(Sol_Min[iVar], field(jPoint, iVar)); - Sol_Max[iVar] = max(Sol_Max[iVar], field(jPoint, iVar)); - } + UpdateMinMax(field[jPoint], rotation); } for (iVar = 0; iVar < ICOUNT; iVar++) { @@ -948,13 +1028,36 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, bufDSend[buf_offset+ICOUNT+iVar] = Sol_Max[iVar]; } - /*--- Rotate the momentum components of the min/max. ---*/ - - if (rotate_periodic) { + /*--- Preserve the original arithmetic for translational periodicity, including its AD recording. ---*/ + if (rotate_periodic && !rotation) { Rotate(zeros, &Sol_Min[1], &bufDSend[buf_offset+1]); Rotate(zeros, &Sol_Max[1], &bufDSend[buf_offset+ICOUNT+1]); } + /*--- The limiters of the velocity components cannot be compared across a rotation + (second phase). Instead, we also send the min and max, over "our" edges, of the + reconstruction increments (computed as in computeLimiters_impl) with the velocity + part rotated, the periodic match then computes its limiters with the complete + stencil. Without rotation these are zero, i.e. they have no effect. ---*/ + + if (commType == PERIODIC_LIM_PRIM_1) { + + for (auto iField = 0u; iField < ICOUNT; iField++) + Sol_Min[iField] = Sol_Max[iField] = 0.0; + + if (rotation) { + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + ReconstructionIncrements(iPoint, jPoint, rotPrim_i); + UpdateMinMax(rotPrim_i, true); + } + } + + for (auto iField = 0u; iField < ICOUNT; iField++) { + bufDSend[buf_offset+2*ICOUNT+iField] = Sol_Min[iField]; + bufDSend[buf_offset+3*ICOUNT+iField] = Sol_Max[iField]; + } + } + break; case PERIODIC_LIM_PRIM_2: @@ -968,7 +1071,14 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, bufDSend[buf_offset+iVar] = limiter(iPoint, iVar); } - if (rotate_periodic) { + /*--- The limiters of the velocity components cannot be compared across a rotation, + they were computed with the complete stencil (see the first phase), we send + a value that is never the minimum. ---*/ + + if (rotation) { + for (auto iCoordinate = 0u; iCoordinate < nDim; iCoordinate++) + bufDSend[buf_offset+1+iCoordinate] = std::numeric_limits::max(); + } else if (rotate_periodic) { Rotate(zeros, &limiter(iPoint,1), &bufDSend[buf_offset+1]); } @@ -1289,6 +1399,19 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, base_nodes->GetSolution_Max()(iPoint, iVar) = Solution_Max; } + /*--- Min/max reconstruction increments over the edges of the periodic match, + used to start the search over "our" edges (only with rotation). ---*/ + + if ((commType == PERIODIC_LIM_PRIM_1) && !PeriodicProj.empty()) { + auto* projections = GetPeriodicProjection(iPoint); + assert(projections != nullptr); + for (auto iField = 0u; iField < ICOUNT; iField++) { + projections[iField] = min(projections[iField], bufDRecv[buf_offset+2*ICOUNT+iField]); + projections[ICOUNT+iField] = max(projections[ICOUNT+iField], + bufDRecv[buf_offset+3*ICOUNT+iField]); + } + } + break; case PERIODIC_LIM_PRIM_2: diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index f63748ba9283..cc24dd1f2601 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -158,7 +158,7 @@ def main(): periodic2d.cfg_dir = "navierstokes/periodic2D" periodic2d.cfg_file = "config.cfg" periodic2d.test_iter = 1400 - periodic2d.test_vals = [-10.817607, -8.363541, -8.287457, -5.334100, -1.088412, -2945.200000] + periodic2d.test_vals = [-10.314545, -8.056213, -7.708669, -4.833054, -1.059806, -2953.400000] test_list.append(periodic2d) ########################## @@ -602,7 +602,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.907561, -12.214137, -19.151426, -13.451780, -19.085218, -13.450796, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.959252, -12.213981, -19.314222, -13.600248, -19.076384, -13.444335, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -627,7 +627,7 @@ def main(): multi_interface.cfg_dir = "turbomachinery/multi_interface" multi_interface.cfg_file = "multi_interface_rst.cfg" multi_interface.test_iter = 5 - multi_interface.test_vals = [-8.632240, -8.894740, -9.348706] + multi_interface.test_vals = [-8.634571, -8.895558, -9.348754] multi_interface.test_vals_aarch64 = [-8.632229, -8.894737, -9.348730] test_list.append(multi_interface) diff --git a/TestCases/navierstokes/periodic2D/multigrid.cfg b/TestCases/navierstokes/periodic2D/multigrid.cfg new file mode 100644 index 000000000000..d7a7c67bbb31 --- /dev/null +++ b/TestCases/navierstokes/periodic2D/multigrid.cfg @@ -0,0 +1,86 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Rotational periodic sector with multigrid, % +% without slope limiter, at lower CFL % +% Author: SU2 Contributors % +% Date: Oct 2026 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +SOLVER= NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 0.1 +AOA= 22.5 +INIT_OPTION= TD_CONDITIONS +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_TEMPERATURE= 300 +FREESTREAM_PRESSURE= 99000 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.00 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1 +REF_AREA= 1 +% +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.87 +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.001 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_PERIODIC= ( per1, per2, 0,0,0, 0,0,45, 0,0,0 ) +MARKER_OUTLET= ( outlet, 99000.0 ) +INLET_TYPE= TOTAL_CONDITIONS +MARKER_INLET= ( inlet, 300, 100000, 1, 0, 0 ) +SPECIFIED_INLET_PROFILE= YES +INLET_FILENAME= inlet.dat +% +MARKER_ANALYZE= ( outlet, inlet ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 20 +CFL_ADAPT= NO +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 0.1 +LINEAR_SOLVER_ITER= 4 + +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 2 +MG_MIN_MESHSIZE= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -11 +CONV_STARTITER= 0 +ITER= 5000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= sector.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_RES, LINSOL_RESIDUAL, SURFACE_PRESSURE_DROP ) diff --git a/TestCases/navierstokes/periodic2D/no_limiter.cfg b/TestCases/navierstokes/periodic2D/no_limiter.cfg new file mode 100644 index 000000000000..2ac3b09dfddb --- /dev/null +++ b/TestCases/navierstokes/periodic2D/no_limiter.cfg @@ -0,0 +1,81 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Rotational periodic sector without slope limiter, % +% testing the periodic Jacobian % +% Author: SU2 Contributors % +% Date: Oct 2026 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +SOLVER= NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 0.1 +AOA= 22.5 +INIT_OPTION= TD_CONDITIONS +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_TEMPERATURE= 300 +FREESTREAM_PRESSURE= 99000 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.00 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1 +REF_AREA= 1 +% +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.87 +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.001 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_PERIODIC= ( per1, per2, 0,0,0, 0,0,45, 0,0,0 ) +MARKER_OUTLET= ( outlet, 99000.0 ) +INLET_TYPE= TOTAL_CONDITIONS +MARKER_INLET= ( inlet, 300, 100000, 1, 0, 0 ) +SPECIFIED_INLET_PROFILE= YES +INLET_FILENAME= inlet.dat +% +MARKER_ANALYZE= ( outlet, inlet ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 100 +CFL_ADAPT= NO +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 0.1 +LINEAR_SOLVER_ITER= 4 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -11 +CONV_STARTITER= 0 +ITER= 5000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= sector.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_RES, LINSOL_RESIDUAL, SURFACE_PRESSURE_DROP ) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index e79219efbeba..ca6e13877235 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -421,6 +421,30 @@ def main(): poiseuille_profile.tol = [0.001, 0.001, 1e-5, 1e-5, 1e-5] test_list.append(poiseuille_profile) + # 2D rotational periodic sector without limiter (Jacobian of the periodic points) + periodic2d_no_limiter = TestCase('periodic2d_no_limiter') + periodic2d_no_limiter.cfg_dir = "navierstokes/periodic2D" + periodic2d_no_limiter.cfg_file = "no_limiter.cfg" + periodic2d_no_limiter.test_iter = 100 + periodic2d_no_limiter.test_vals = [-3.216881, -0.582813, -0.629049, 2.266141, -1.056877, -811.470000] + test_list.append(periodic2d_no_limiter) + + # 2D rotational periodic sector with multigrid + periodic2d_multigrid = TestCase('periodic2d_multigrid') + periodic2d_multigrid.cfg_dir = "navierstokes/periodic2D" + periodic2d_multigrid.cfg_file = "multigrid.cfg" + periodic2d_multigrid.test_iter = 300 + periodic2d_multigrid.test_vals = [-4.422778, -1.615251, -1.453736, 1.047350, -1.410225, -2060.600000] + test_list.append(periodic2d_multigrid) + + # 2D rotational periodic sector with limiter + periodic2d = TestCase('periodic2d') + periodic2d.cfg_dir = "navierstokes/periodic2D" + periodic2d.cfg_file = "config.cfg" + periodic2d.test_iter = 100 + periodic2d.test_vals = [-3.266248, -0.617784, -0.620252, 2.216476, -1.007761, -1035.000000] + test_list.append(periodic2d) + ########################## ### Compressible RANS ### ########################## @@ -1252,7 +1276,7 @@ def main(): Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 Aachen_3D_restart.tol = 0.00001 - Aachen_3D_restart.test_vals = [-7.701420, -8.504728, -6.014939, -6.468223, -5.801124, -4.607179, -5.550665, -5.300778, -3.804188, -5.255983, -5.763060, -3.609605, -2.229249, -2.880453, -0.563469] + Aachen_3D_restart.test_vals = [-7.688483, -8.466346, -6.035067, -6.425967, -5.820125, -4.597817, -5.528297, -5.301466, -3.819862, -5.242356, -5.738805, -3.630202, -2.213240, -2.857919, -0.576864] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart @@ -1260,7 +1284,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.941917, -12.212515, -19.254664, -13.545311, -19.087161, -13.454459, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.905382, -12.212100, -19.212143, -13.508268, -19.078783, -13.446389, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -1285,7 +1309,7 @@ def main(): multi_interface.cfg_dir = "turbomachinery/multi_interface" multi_interface.cfg_file = "multi_interface_rst.cfg" multi_interface.test_iter = 5 - multi_interface.test_vals = [-8.632227, -8.894736, -9.348706] + multi_interface.test_vals = [-8.634558, -8.895554, -9.348754] multi_interface.test_vals_aarch64 = [-8.632227, -8.894736, -9.348706] test_list.append(multi_interface) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 26a32c3e2762..b80601e77efe 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -251,7 +251,7 @@ def main(): discadj_trans_stator.cfg_dir = "disc_adj_turbomachinery/transonic_stator_2D" discadj_trans_stator.cfg_file = "transonic_stator.cfg" discadj_trans_stator.test_iter = 79 - discadj_trans_stator.test_vals = [79.000000, -7.555663, -10.335501, -10.356934, -13.629559] + discadj_trans_stator.test_vals = [79.000000, -7.555647, -10.335486, -10.356919, -13.629543] discadj_trans_stator.test_vals_aarch64 = [79.000000, -7.555647, -10.335486, -10.356919, -13.629543] test_list.append(discadj_trans_stator) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index e1a826cbf170..f227b308e734 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -949,7 +949,7 @@ def main(): Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 - Aachen_3D_restart.test_vals = [-7.701421, -8.504727, -6.014939, -6.468221, -5.801125, -4.607173, -5.550665, -5.300779, -3.804187, -5.255982, -5.763060, -3.609601, -2.229250, -2.880453, -0.563470] + Aachen_3D_restart.test_vals = [-7.689117, -8.476204, -6.036107, -6.427402, -5.831576, -4.599024, -5.530512, -5.305687, -3.821789, -5.243255, -5.754831, -3.631603, -2.213934, -2.895728, -0.578545] Aachen_3D_restart.enabled_with_asan = False test_list.append(Aachen_3D_restart) @@ -958,7 +958,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.944235, -12.212620, -19.261137, -13.549357, -19.083828, -13.444697, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-11.905642, -12.212181, -19.212554, -13.508053, -19.069991, -13.439494, 73286.000000, 73286.000000, 0.020056, 82.286000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -983,7 +983,7 @@ def main(): multi_interface.cfg_dir = "turbomachinery/multi_interface" multi_interface.cfg_file = "multi_interface_rst.cfg" multi_interface.test_iter = 5 - multi_interface.test_vals = [-8.632227, -8.894736, -9.348706] + multi_interface.test_vals = [-8.634558, -8.895554, -9.348754] multi_interface.test_vals_aarch64 = [-8.632227, -8.894736, -9.348706] test_list.append(multi_interface) diff --git a/TestCases/turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg b/TestCases/turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg index 0aebb8cd4d42..9e1c0f04cc45 100755 --- a/TestCases/turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg +++ b/TestCases/turbomachinery/Aachen_turbine/aachen_3D_MP_restart.cfg @@ -228,6 +228,10 @@ LINEAR_SOLVER_ERROR= 1E-4 % Max number of iterations of the linear solver for the implicit formulation LINEAR_SOLVER_ITER= 15 % +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 1 +% % ----------------------- SLOPE LIMITER DEFINITION ----------------------------% % % Coefficient for the limiter diff --git a/UnitTests/SU2_CFD/periodic_limiters.cpp b/UnitTests/SU2_CFD/periodic_limiters.cpp new file mode 100644 index 000000000000..22b4026e9374 --- /dev/null +++ b/UnitTests/SU2_CFD/periodic_limiters.cpp @@ -0,0 +1,148 @@ +/*! + * \file periodic_limiters.cpp + * \brief Tests for periodic limiter storage and reconstruction helpers. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include "../../SU2_CFD/include/solvers/CSolver.hpp" + +TEST_CASE("Configured zero angles select translation and tiny angles select rotation", "[PeriodicLimiter]") { + const su2double noAngles[3] = {0.0, -0.0, 0.0}, negativeAngle[3] = {0.0, -0.5, 0.0}; + CHECK_FALSE(GeometryToolbox::HasRotation(noAngles)); + CHECK(GeometryToolbox::HasRotation(negativeAngle)); + const su2double deg2rad = PI_NUMBER / 180.0; + for (int iDim = 0; iDim < 3; ++iDim) { + su2double angles[3] = {0.0, 0.0, 0.0}; + angles[iDim] = su2double(0.0) * deg2rad; + CHECK_FALSE(GeometryToolbox::HasRotation(angles)); + angles[iDim] *= -1.0; + CHECK_FALSE(GeometryToolbox::HasRotation(angles)); + angles[iDim] = su2double(1.0e-18) * deg2rad; + CHECK(GeometryToolbox::HasRotation(angles)); + angles[iDim] *= -1.0; + CHECK(GeometryToolbox::HasRotation(angles)); + } +} + +TEST_CASE("Rotation of component bounds encloses every corner", "[PeriodicLimiter]") { + auto checkCorners = [](auto& rotation, auto& lower, auto& upper) { + constexpr int nDim = sizeof(lower) / sizeof(lower[0]); + su2double expectedMin[nDim], expectedMax[nDim]; + for (int iDim = 0; iDim < nDim; ++iDim) { + expectedMin[iDim] = std::numeric_limits::max(); + expectedMax[iDim] = -expectedMin[iDim]; + } + const su2double origin[nDim] = {0.0}; + for (int corner = 0; corner < (1 << nDim); ++corner) { + su2double point[nDim], rotated[nDim]; + for (int iDim = 0; iDim < nDim; ++iDim) point[iDim] = (corner & (1 << iDim)) ? upper[iDim] : lower[iDim]; + GeometryToolbox::Rotate(rotation, origin, point, rotated); + for (int iDim = 0; iDim < nDim; ++iDim) { + expectedMin[iDim] = std::min(expectedMin[iDim], rotated[iDim]); + expectedMax[iDim] = std::max(expectedMax[iDim], rotated[iDim]); + } + } + GeometryToolbox::RotateBox(rotation, lower, upper); + for (int iDim = 0; iDim < nDim; ++iDim) { + CHECK(lower[iDim] == Approx(expectedMin[iDim])); + CHECK(upper[iDim] == Approx(expectedMax[iDim])); + } + }; + + SECTION("2D oblique rotation") { + su2double rotation[2][2], lower[2] = {-2.0, 1.0}, upper[2] = {3.0, 4.0}; + GeometryToolbox::RotationMatrix(su2double(PI_NUMBER / 4.0), rotation); + checkCorners(rotation, lower, upper); + } + SECTION("3D rotation about all axes") { + su2double rotation[3][3], lower[3] = {-2.0, 1.0, -3.0}, upper[3] = {3.0, 4.0, 2.0}; + GeometryToolbox::RotationMatrix(su2double(0.3), su2double(-0.7), su2double(1.2), rotation); + checkCorners(rotation, lower, upper); + } +} + +TEST_CASE("Shared reconstruction increment has the MUSCL scaling", "[PeriodicLimiter]") { + const su2double halfEdge[3] = {1.0, -1.0, 0.5}; + const su2double gradient[3] = {2.0, -1.0, 3.0}; + for (const auto kappa : {-1.0, 0.0, 0.5, 1.0}) { + const auto increment = LimiterHelpers<>::reconstructionIncrement(3, halfEdge, gradient, su2double(7.0), + su2double(15.0), su2double(kappa)); + CHECK(increment == Approx(4.5 - 0.5 * kappa)); + CHECK(LimiterHelpers<>::reconstructionIncrement(2, halfEdge, gradient, su2double(7.0), su2double(15.0), + su2double(kappa)) == Approx(3.0 + kappa)); + CHECK(LimiterHelpers<>::reconstructionIncrement(2, halfEdge, gradient, su2double(7.0), su2double(13.0), + su2double(kappa)) == Approx(3.0)); + const auto linearIncrement = LimiterHelpers<>::reconstructionIncrement(3, halfEdge, gradient, su2double(7.0), + su2double(16.0), su2double(kappa)); + CHECK(linearIncrement == Approx(4.5)); + } +} + +TEST_CASE("Periodic projection storage uses unique receive points", "[PeriodicLimiter]") { + std::stringstream options; + options << "SOLVER= EULER\n" + "MARKER_PERIODIC= (per1, per2, 0,0,0, 0,0,45, 0,0,0)\n"; + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + config.SetnMarker_All(1); + config.SetMarker_All_KindBC(0, PERIODIC_BOUNDARY); + config.SetMarker_All_TagBound(0, "per1"); + + struct ProjectionSolver : CSolver { + CVariable* GetBaseClassPointerToNodes() override { return nullptr; } + ProjectionSolver() { + nPoint = 1000000; + nVar = 4; + nPrimVarGrad = 4; + rotate_periodic = true; + } + } solver; + + CGeometry geometry; + geometry.nPeriodicRecv = 1; + geometry.nPoint_PeriodicRecv = new int[2]{0, 4}; + geometry.Local_Point_PeriodicRecv = new unsigned long[4]{12, 27, 12, 99999}; + + auto* storage = solver.GetPeriodicProjections(geometry, config); + REQUIRE(storage != nullptr); + CHECK(storage->rows() == 3); + CHECK(storage->cols() == 8); + CHECK(solver.GetPeriodicProjection(13) == nullptr); + REQUIRE(solver.GetPeriodicProjection(12) != nullptr); + REQUIRE(solver.GetPeriodicProjection(27) != nullptr); + REQUIRE(solver.GetPeriodicProjection(99999) != nullptr); + solver.GetPeriodicProjection(12)[0] = -2.0; + CHECK(solver.GetPeriodicProjection(27) != solver.GetPeriodicProjection(12)); + CHECK(solver.GetPeriodicProjections(geometry, config) == storage); + CHECK(solver.GetPeriodicProjection(12)[0] == -2.0); + + CGeometry noReceives; + ProjectionSolver solverWithoutReceives; + auto* emptyStorage = solverWithoutReceives.GetPeriodicProjections(noReceives, config); + REQUIRE(emptyStorage != nullptr); + CHECK(emptyStorage->empty()); + CHECK(solverWithoutReceives.GetPeriodicProjection(12) == nullptr); + + solver.SetRotatePeriodic(false); + CHECK(solver.GetPeriodicProjections(geometry, config) == nullptr); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index ba3c63afce92..1d1a97f3b827 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -19,6 +19,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'SU2_CFD/output/CNEMOCompOutput_tests.cpp', 'SU2_CFD/output/COutput_convergence_tests.cpp', 'SU2_CFD/gradients.cpp', + 'SU2_CFD/periodic_limiters.cpp', 'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp', 'SU2_CFD/nemo_viscous_assembly.cpp', 'SU2_CFD/windowing.cpp',