From e4deaa2ee5497ef3f84e64cb8863151d7765bb89 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:33 +0200 Subject: [PATCH 01/10] Stop when streamwise periodicity is used with multigrid The run crashed before the first iteration, the streamwise periodic quantities are only set up on the fine grid. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- Common/src/CConfig.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index fd4207b0c5f..b9a9d8eea95 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -5557,6 +5557,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("Streamwise Periodicity only works with \"INC_NONDIM= DIMENSIONAL\", the nondimensionalization with source terms doesn;t work in general.", CURRENT_FUNCTION); if (Axisymmetric) SU2_MPI::Error("Streamwise Periodicity terms does not not have axisymmetric corrections.", CURRENT_FUNCTION); + if (nMGLevels != 0) + SU2_MPI::Error("Streamwise Periodicity does not support multigrid, use MGLEVEL= 0.", CURRENT_FUNCTION); if (!Energy_Equation) Streamwise_Periodic_Temperature = false; } else { /*--- Safety measure ---*/ From 071ad4d045ab3c0cbabc863b3e6d126c2bcecdbe Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:34 +0200 Subject: [PATCH 02/10] Use INTEGRATED_HEATFLUX in the integrated heat of streamwise periodic flow With INTEGRATED_HEATFLUX= YES the value of MARKER_HEATFLUX is in W, it was multiplied by the area again as if it was in W/m2. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- SU2_CFD/src/solvers/CIncNSSolver.cpp | 3 +- .../sp_pipeSlice_3d_dp_ihf_tp.cfg | 85 +++++++++++++++++++ TestCases/parallel_regression.py | 8 ++ 3 files changed, 95 insertions(+), 1 deletion(-) create mode 100644 TestCases/incomp_navierstokes/streamwise_periodic/pipeSlice_3d/sp_pipeSlice_3d_dp_ihf_tp.cfg diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index ec1b78935e1..7148b2d3715 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -209,7 +209,8 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry, /*--- Identify the boundary by string name and retrive heatflux from config ---*/ const auto Marker_StringTag = config->GetMarker_All_TagBound(iMarker); - const su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_StringTag); + su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_StringTag); + if (config->GetIntegrated_HeatFlux()) Wall_HeatFlux /= geometry->GetSurfaceArea(config, iMarker); for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) { diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/pipeSlice_3d/sp_pipeSlice_3d_dp_ihf_tp.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/pipeSlice_3d/sp_pipeSlice_3d_dp_ihf_tp.cfg new file mode 100644 index 00000000000..37d1ff9a70a --- /dev/null +++ b/TestCases/incomp_navierstokes/streamwise_periodic/pipeSlice_3d/sp_pipeSlice_3d_dp_ihf_tp.cfg @@ -0,0 +1,85 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Same as sp_pipeSlice_3d_dp_hf_tp.cfg with the energy % +% equation, periodic temperature and a wall heat given as % +% integrated value (INTEGRATED_HEATFLUX) % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_NAVIER_STOKES +% +KIND_TURB_MODEL= NONE +% +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_ENERGY_EQUATION = YES +INC_DENSITY_INIT= 1.0 +INC_VELOCITY_INIT= ( 0.0, 0.0, 0.3 ) +INC_NONDIM= DIMENSIONAL +% +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.8e-5 +% +% --------------------- STREAMWISE PERIODICITY DEFINITION ---------------------% +% +KIND_STREAMWISE_PERIODIC= PRESSURE_DROP +STREAMWISE_PERIODIC_PRESSURE_DROP= 0.001 +STREAMWISE_PERIODIC_TEMPERATURE= YES +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +% Heat [W] over the whole wall (area 1.570028349e-5 m2, i.e. 1000 W/m2) +MARKER_HEATFLUX= (wall, 0.01570028349) +INTEGRATED_HEATFLUX= YES +MARKER_PERIODIC= ( inlet, outlet, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0005 ) +% +MARKER_PLOTTING= ( inlet ) +MARKER_MONITORING= ( wall ) +MARKER_ANALYZE = ( oulet ) +MARKER_ANALYZE_AVERAGE = AREA +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +%CFL_NUMBER= 1e10 +CFL_NUMBER= 50000 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 0.5, 10, 15.0, 1e30 ) +ITER= 15000 +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-15 +LINEAR_SOLVER_ITER= 10 +% +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -24 +CONV_STARTITER= 10 +% +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= pipe1cell3D.su2 +% +OUTPUT_FILES= ( RESTART, PARAVIEW_MULTIBLOCK, SURFACE_TECPLOT_ASCII ) +OUTPUT_WRT_FREQ= 1000 +% +HISTORY_OUTPUT= ( RMS_RES, FLOW_COEFF, STREAMWISE_PERIODIC, LINSOL ) +SCREEN_OUTPUT= ( INNER_ITER, RMS_PRESSURE, RMS_ENTHALPY, STREAMWISE_MASSFLOW, STREAMWISE_HEAT ) +SCREEN_WRT_FREQ_INNER= 100 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index e85b5b7be66..71f21c36575 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -772,6 +772,14 @@ def main(): sp_pipeSlice_3d_dp_hf_tp.test_vals = [-11.119796, -11.234737, -9.429479, -0.000023] test_list.append(sp_pipeSlice_3d_dp_hf_tp) + # Same with the energy equation, periodic temperature and an integrated wall heat + sp_pipeSlice_3d_dp_ihf_tp = TestCase('sp_pipeSlice_3d_dp_ihf_tp') + sp_pipeSlice_3d_dp_ihf_tp.cfg_dir = "incomp_navierstokes/streamwise_periodic/pipeSlice_3d" + sp_pipeSlice_3d_dp_ihf_tp.cfg_file = "sp_pipeSlice_3d_dp_ihf_tp.cfg" + sp_pipeSlice_3d_dp_ihf_tp.test_iter = 10 + sp_pipeSlice_3d_dp_ihf_tp.test_vals = [-11.120123, -5.196299, -0.000023, -0.015700] + test_list.append(sp_pipeSlice_3d_dp_ihf_tp) + # 2D pin array with heat transfer BC on pin surfaces inc_heatTransfer_BC = TestCase('inc_heatTransfer_BC') inc_heatTransfer_BC.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d" From 898254959e71cd108df4367d6bf157061beb3a18 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:34 +0200 Subject: [PATCH 03/10] Fix the turbulent source term of the streamwise periodic temperature The term with the gradient of the eddy viscosity was divided by the length of the translation vector instead of its square. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index e3797c7c5a5..02270d9f414 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -794,8 +794,10 @@ CNumerics::ResidualType<> CSourceIncStreamwise_Periodic::ComputeResidual(const C /*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/ if(turbulent) { - /*--- Compute a scalar factor ---*/ - scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb); + /*--- Compute a scalar factor. The temperature is T = T_periodic + Q / (mdot * cp * |t|^2) * dot_prod(t*x) + (Patankar et al., J. Heat Transfer 99(2), 1977), so the turbulent heat flux cp * mu_t / Pr_t * grad(T) + has the extra divergence Q / (mdot * |t|^2 * Pr_t) * dot_prod(t*grad(mu_t)). ---*/ + scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation * Prandtl_Turb); /*--- Compute scalar product between periodic translation vector and eddy viscosity gradient. ---*/ dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, AuxVar_Grad_i[0]); From a7a424b2d6d2de6222efa2b7c965ba03f7d75672 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:34 +0200 Subject: [PATCH 04/10] Use the signed distance to recover the streamwise periodic pressure and temperature The sum of the absolute values of each component is only the distance along the translation vector if the vector is aligned with an axis. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- SU2_CFD/src/solvers/CIncNSSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 7148b2d3715..eb4ef77dea9 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -256,7 +256,7 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config, /*--- First, compute helping terms based on relative distance (0,l) between periodic markers ---*/ su2double dot_product = 0.0; for (unsigned short iDim = 0; iDim < nDim; iDim++) - dot_product += fabs( (geometry->nodes->GetCoord(iPoint,iDim) - ReferenceNode[iDim]) * config->GetPeriodic_Translation(0)[iDim]); + dot_product += (geometry->nodes->GetCoord(iPoint,iDim) - ReferenceNode[iDim]) * config->GetPeriodic_Translation(0)[iDim]; /*--- Second, substract/add correction from reduced pressure/temperature to get recoverd pressure/temperature ---*/ const su2double Pressure_Recovered = nodes->GetPressure(iPoint) - SPvals.Streamwise_Periodic_PressureDrop / From 192b32783c7f740e9960d3458a84a09b83cba8dc Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 6 Oct 2026 01:00:34 +0200 Subject: [PATCH 05/10] Initialize the streamwise periodic values The heat was printed with whatever was in memory when the energy equation was off. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- Common/include/option_structure.hpp | 12 ++++++------ SU2_CFD/src/solvers/CIncNSSolver.cpp | 2 +- TestCases/parallel_regression.py | 4 ++-- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bcd2e9c0632..02c079b8cac 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2989,12 +2989,12 @@ static const MapType Streamwise_Periodic_ * \brief Container to hold Variables for streamwise Periodic flow as they are often used together in places. */ struct StreamwisePeriodicValues { - su2double Streamwise_Periodic_PressureDrop; /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */ - su2double Streamwise_Periodic_MassFlow; /*!< \brief Value of current massflow [kg/s] which results in a delta p and therefore an artificial body force vector. */ - su2double Streamwise_Periodic_IntegratedHeatFlow; /*!< \brief Value of of the net sum of heatflow [W] into the domain. */ - su2double Streamwise_Periodic_InletTemperature; /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */ - su2double Streamwise_Periodic_BoundaryArea; /*!< \brief Global Surface area of the streamwise periodic interface. */ - su2double Streamwise_Periodic_AvgDensity; /*!< \brief Area avg density on the periodic interface. */ + su2double Streamwise_Periodic_PressureDrop = 0.0; /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */ + su2double Streamwise_Periodic_MassFlow = 0.0; /*!< \brief Value of current massflow [kg/s] which results in a delta p and therefore an artificial body force vector. */ + su2double Streamwise_Periodic_IntegratedHeatFlow = 0.0; /*!< \brief Value of of the net sum of heatflow [W] into the domain. */ + su2double Streamwise_Periodic_InletTemperature = 0.0; /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */ + su2double Streamwise_Periodic_BoundaryArea = 0.0; /*!< \brief Global Surface area of the streamwise periodic interface. */ + su2double Streamwise_Periodic_AvgDensity = 0.0; /*!< \brief Area avg density on the periodic interface. */ }; /*! diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index eb4ef77dea9..0042b62b0f2 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -255,7 +255,7 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config, /*--- First, compute helping terms based on relative distance (0,l) between periodic markers ---*/ su2double dot_product = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) + for (auto iDim = 0u; iDim < nDim; iDim++) dot_product += (geometry->nodes->GetCoord(iPoint,iDim) - ReferenceNode[iDim]) * config->GetPeriodic_Translation(0)[iDim]; /*--- Second, substract/add correction from reduced pressure/temperature to get recoverd pressure/temperature ---*/ diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 71f21c36575..6aded1ee082 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -774,8 +774,8 @@ def main(): # Same with the energy equation, periodic temperature and an integrated wall heat sp_pipeSlice_3d_dp_ihf_tp = TestCase('sp_pipeSlice_3d_dp_ihf_tp') - sp_pipeSlice_3d_dp_ihf_tp.cfg_dir = "incomp_navierstokes/streamwise_periodic/pipeSlice_3d" - sp_pipeSlice_3d_dp_ihf_tp.cfg_file = "sp_pipeSlice_3d_dp_ihf_tp.cfg" + sp_pipeSlice_3d_dp_ihf_tp.cfg_dir = "incomp_navierstokes/streamwise_periodic/pipeSlice_3d" + sp_pipeSlice_3d_dp_ihf_tp.cfg_file = "sp_pipeSlice_3d_dp_ihf_tp.cfg" sp_pipeSlice_3d_dp_ihf_tp.test_iter = 10 sp_pipeSlice_3d_dp_ihf_tp.test_vals = [-11.120123, -5.196299, -0.000023, -0.015700] test_list.append(sp_pipeSlice_3d_dp_ihf_tp) From 0c8e453ecc70d86fa18e95b89460be6581c3d883 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 10:31:33 +0200 Subject: [PATCH 06/10] Update streamwise periodic pressure drop once per OpenMP team --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 30 +++++++++------- TestCases/hybrid_regression.py | 8 +++++ .../streamwise_periodic/massflow_box.cfg | 36 +++++++++++++++++++ 3 files changed, 61 insertions(+), 13 deletions(-) create mode 100644 TestCases/incomp_navierstokes/streamwise_periodic/massflow_box.cfg diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 512a99a97dc..a2bbef972cb 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -1889,19 +1889,23 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont /*--- The Pressure drop is iteratively adapted to result in the prescribed Target-Massflow. ---*/ /*---------------------------------------------------------------------------------------------*/ - /*--- Compute update to Delta p based on massflow-difference ---*/ - const su2double Average_Density_Global = SPvals.Streamwise_Periodic_AvgDensity; - const su2double Area_Global = SPvals.Streamwise_Periodic_BoundaryArea; - const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref()); - const su2double MassFlow_Global = SPvals.Streamwise_Periodic_MassFlow; - const su2double ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2)); - - /*--- Store updated pressure difference ---*/ - const su2double damping_factor = config->GetInc_Outlet_Damping(); - SPvalsUpdated = SPvals; - SPvalsUpdated.Streamwise_Periodic_PressureDrop += damping_factor*ddP; - if (!config->GetDiscrete_Adjoint()) - SPvals = SPvalsUpdated; + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + { + /*--- Compute update to Delta p based on massflow-difference ---*/ + const su2double Average_Density_Global = SPvals.Streamwise_Periodic_AvgDensity; + const su2double Area_Global = SPvals.Streamwise_Periodic_BoundaryArea; + const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref()); + const su2double MassFlow_Global = SPvals.Streamwise_Periodic_MassFlow; + const su2double ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2)); + + /*--- Store updated pressure difference ---*/ + const su2double damping_factor = config->GetInc_Outlet_Damping(); + SPvalsUpdated = SPvals; + SPvalsUpdated.Streamwise_Periodic_PressureDrop += damping_factor*ddP; + if (!config->GetDiscrete_Adjoint()) + SPvals = SPvalsUpdated; + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS /*--- Set delta_p, m_dot, inlet_T, integrated_heat ---*/ numerics->SetStreamwisePeriodicValues(SPvalsUpdated); diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index f63748ba928..c7cbd6a9753 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -153,6 +153,14 @@ def main(): poiseuille_profile.test_vals_aarch64 = [-12.004276, -7.636719, -0.000000, 2.089953] test_list.append(poiseuille_profile) + # Streamwise periodic pressure-drop update with two OpenMP threads. + sp_massflow_box = TestCase('sp_massflow_box') + sp_massflow_box.cfg_dir = "incomp_navierstokes/streamwise_periodic" + sp_massflow_box.cfg_file = "massflow_box.cfg" + sp_massflow_box.test_iter = 4 + sp_massflow_box.test_vals = [-4.914357, 0.248162, -0.092770] + test_list.append(sp_massflow_box) + # 2D Rotational Periodic periodic2d = TestCase('periodic2d') periodic2d.cfg_dir = "navierstokes/periodic2D" diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/massflow_box.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/massflow_box.cfg new file mode 100644 index 00000000000..8763d36ff6d --- /dev/null +++ b/TestCases/incomp_navierstokes/streamwise_periodic/massflow_box.cfg @@ -0,0 +1,36 @@ +% Small generated periodic channel: pressure-drop update must be independent of OpenMP threads. +SOLVER= INC_NAVIER_STOKES +MESH_FORMAT= BOX +MESH_BOX_SIZE= (9, 9, 9) +MESH_BOX_LENGTH= (1, 1, 1) +MESH_BOX_OFFSET= (0, 0, 0) +INC_DENSITY_MODEL= CONSTANT +INC_DENSITY_INIT= 1.0 +INC_VELOCITY_INIT= (0.1, 0.0, 0.0) +INC_TEMPERATURE_INIT= 300.0 +INC_ENERGY_EQUATION= NO +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.01 +MARKER_HEATFLUX= (y_minus,0, y_plus,0, z_minus,0, z_plus,0) +MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0) +KIND_STREAMWISE_PERIODIC= MASSFLOW +STREAMWISE_PERIODIC_MASSFLOW= 1.0 +STREAMWISE_PERIODIC_PRESSURE_DROP= 0.0 +INC_OUTLET_DAMPING= 0.1 +NUM_METHOD_GRAD= GREEN_GAUSS +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= NO +TIME_DISCRE_FLOW= EULER_IMPLICIT +CFL_NUMBER= 1.0 +CFL_ADAPT= NO +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ITER= 10 +ITER= 5 +CONV_RESIDUAL_MINVAL= -30 +OUTPUT_FILES= NONE +OUTPUT_WRT_FREQ= 5 +HISTORY_OUTPUT= (RMS_RES, STREAMWISE_PERIODIC) +SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, STREAMWISE_DP, STREAMWISE_MASSFLOW) +SCREEN_WRT_FREQ_INNER= 1 +INC_NONDIM= DIMENSIONAL From 609d20757af2a014fa6e5208af70db5f5591c8d5 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 10:31:33 +0200 Subject: [PATCH 07/10] Complete periodic Green-Gauss auxiliary gradients --- Common/include/option_structure.hpp | 1 + SU2_CFD/src/solvers/CSolver.cpp | 18 +++++++++++--- UnitTests/SU2_CFD/gradients.cpp | 38 +++++++++++++++++++++++++++++ 3 files changed, 54 insertions(+), 3 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 02c079b8cac..f75e76aae66 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2853,6 +2853,7 @@ enum PERIODIC_QUANTITIES { PERIODIC_LIM_PRIM_1 , /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */ PERIODIC_LIM_PRIM_2 , /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */ PERIODIC_IMPLICIT , /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */ + PERIODIC_AUXVAR_GG , /*!< \brief Auxiliary scalar gradient communication for Green-Gauss. */ }; /*! diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 56cee4b6ad5..b6897ef98e1 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -229,6 +229,12 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, COUNT_PER_POINT = 2; MPI_TYPE = COMM_TYPE::DOUBLE; break; + case PERIODIC_AUXVAR_GG: + COUNT_PER_POINT = base_nodes->GetnAuxVar()*nDim; + MPI_TYPE = COMM_TYPE::DOUBLE; + ICOUNT = base_nodes->GetnAuxVar(); + JCOUNT = nDim; + break; case PERIODIC_SOL_GG: case PERIODIC_SOL_GG_R: COUNT_PER_POINT = nVar*nDim; @@ -296,6 +302,8 @@ namespace PeriodicCommHelpers { case PERIODIC_PRIM_ULS: return nodes->GetGradient_Primitive(); break; + case PERIODIC_AUXVAR_GG: + return nodes->GetAuxVarGradient(); case PERIODIC_SOL_GG: case PERIODIC_SOL_LS: case PERIODIC_SOL_ULS: @@ -309,6 +317,8 @@ namespace PeriodicCommHelpers { const su2activematrix& selectField(CVariable* nodes, unsigned short commType) { switch(commType) { + case PERIODIC_AUXVAR_GG: + return nodes->GetAuxVar(); case PERIODIC_PRIM_GG: case PERIODIC_PRIM_LS: case PERIODIC_PRIM_ULS: @@ -725,6 +735,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, case PERIODIC_SOL_GG_R: case PERIODIC_PRIM_GG: case PERIODIC_PRIM_GG_R: + case PERIODIC_AUXVAR_GG: /*--- Access and rotate the partial G-G gradient. These will be summed on both sides of the periodic faces before dividing @@ -742,9 +753,9 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, Rotate(zeros, jacBlock[iVar], rotBlock[iVar]); } - /*--- Rotate the vector components of the solution. ---*/ + /*--- Auxiliary fields are scalars: rotate their spatial gradients only. ---*/ - if (rotate_periodic) { + if (rotate_periodic && commType != PERIODIC_AUXVAR_GG) { for (iDim = 0; iDim < nDim; iDim++) { su2double d_diDim[3] = {0.0}; for (iVar = 1; iVar < 1+nDim; ++iVar) { @@ -1233,6 +1244,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, case PERIODIC_SOL_GG_R: case PERIODIC_PRIM_GG: case PERIODIC_PRIM_GG_R: + case PERIODIC_AUXVAR_GG: /*--- For G-G, we accumulate partial gradients then compute the final value using the entire volume of the periodic cell. ---*/ @@ -2204,7 +2216,7 @@ void CSolver::SetAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config) const auto& solution = base_nodes->GetAuxVar(); auto& gradient = base_nodes->GetAuxVarGradient(); - computeGradientsGreenGauss(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, + computeGradientsGreenGauss(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_AUXVAR_GG, *geometry, *config, solution, 0, base_nodes->GetnAuxVar(), -1, gradient); } diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index 72baf525ca6..882e363ce5b 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -31,6 +31,7 @@ #include "../../SU2_CFD/include/solvers/CSolver.hpp" #include "../../SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp" #include "../../SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp" +#include "../UnitQuadTestCase.hpp" /*! * \brief Base class for gradient tests using a unit cube geometry. @@ -154,3 +155,40 @@ TEST_CASE("GG", "[Gradients]") { testGreenGauss(); } TEST_CASE("LS", "[Gradients]") { testLeastSquares(false); } TEST_CASE("WLS", "[Gradients]") { testLeastSquares(true); } + +TEST_CASE("Periodic auxiliary Green-Gauss gradient", "[Gradients][Periodic]") { + const bool rotation = GENERATE(false, true); + UnitQuadTestCase field; + const string custom = "MARKER_CUSTOM= ( x_minus, x_plus, z_plus, z_minus)"; + field.config_options.replace(field.config_options.find(custom), custom.size(), "MARKER_CUSTOM= (z_plus, z_minus)"); + field.AddOption(rotation ? "MARKER_PERIODIC= (x_minus,x_plus, 0,0.5,0.5, 90,0,0, 1,0,0)" + : "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); + field.AddOption("VORTICITY_CONFINEMENT= YES"); + field.InitConfig(); + field.InitGeometry(); + 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(); + REQUIRE(nodes->GetnAuxVar() == 1); + for (auto iPoint = 0ul; iPoint < field.geometry->GetnPoint(); ++iPoint) { + const auto* coord = field.geometry->nodes->GetCoord(iPoint); + const auto y = coord[1] - 0.5; + const auto z = coord[2] - 0.5; + nodes->SetAuxVar(iPoint, 0, rotation ? 2.0 + y * y + z * z : 2.0 + 3.0 * coord[1] - coord[2]); + } + solver->SetAuxVar_Gradient_GG(field.geometry.get(), field.config.get()); + + su2double error = 0.0; + for (auto iPoint = 0ul; iPoint < field.geometry->GetnPointDomain(); ++iPoint) { + const auto* coord = field.geometry->nodes->GetCoord(iPoint); + if (rotation && (coord[1] == 0.0 || coord[1] == 1.0 || coord[2] == 0.0 || coord[2] == 1.0)) continue; + const su2double expected[3] = {0.0, rotation ? 2.0 * (coord[1] - 0.5) : 3.0, + rotation ? 2.0 * (coord[2] - 0.5) : -1.0}; + for (auto iDim = 0u; iDim < 3; ++iDim) + error = max(error, abs(nodes->GetAuxVarGradient(iPoint, 0, iDim) - expected[iDim])); + } + CHECK(error < 1e-9); +} From cbf1003cf955a1b4d706204ad3c177ee9a661d89 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:34:21 +0200 Subject: [PATCH 08/10] Exchange auxiliary periodic gradients and correct streamwise heated-wall support --- Common/include/option_structure.hpp | 1 + Common/src/CConfig.cpp | 15 ++ SU2_CFD/src/solvers/CIncNSSolver.cpp | 3 +- SU2_CFD/src/solvers/CSolver.cpp | 24 ++- .../variables/CSpeciesFlameletVariable.cpp | 1 + UnitTests/Common/CConfig_tests.cpp | 45 +++++ UnitTests/Common/check_periodic_support.py | 50 +++++ UnitTests/SU2_CFD/gradients.cpp | 183 +++++++++++++++++- UnitTests/UnitQuadTestCase.hpp | 17 +- UnitTests/meson.build | 3 + 10 files changed, 328 insertions(+), 14 deletions(-) create mode 100644 UnitTests/Common/check_periodic_support.py diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index f75e76aae66..aa124f59a20 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2853,6 +2853,7 @@ enum PERIODIC_QUANTITIES { PERIODIC_LIM_PRIM_1 , /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */ PERIODIC_LIM_PRIM_2 , /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */ PERIODIC_IMPLICIT , /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */ + PERIODIC_AUXVAR_LS , /*!< \brief Auxiliary scalar gradient communication for weighted least-squares. */ PERIODIC_AUXVAR_GG , /*!< \brief Auxiliary scalar gradient communication for Green-Gauss. */ }; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index b9a9d8eea95..b323a51500d 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -37,6 +37,7 @@ #include "../include/basic_types/ad_structure.hpp" #include "../include/toolboxes/printing_toolbox.hpp" +#include "../include/toolboxes/geometry_toolbox.hpp" #include "../include/toolboxes/SwapBytes.hpp" using namespace PrintingToolbox; @@ -5551,6 +5552,20 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("Streamwise Periodic Flow + Incompressible Euler: Not tested yet.", CURRENT_FUNCTION); if (nMarker_PerBound == 0) SU2_MPI::Error("A MARKER_PERIODIC pair has to be set with KIND_STREAMWISE_PERIODIC != NONE.", CURRENT_FUNCTION); + /*--- The first pair defines the streamwise direction used by every source. ---*/ + const auto* translation = GetPeriodic_Translation(0); + if (GeometryToolbox::SquaredNorm(3, translation) == 0) + SU2_MPI::Error("The first streamwise periodic pair must have a nonzero translation.", CURRENT_FUNCTION); + const auto* angles = GetPeriodicRotAngles(Marker_PerBound[0]); + su2double rotation[3][3]; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation); + for (auto i = 0u; i < 3; ++i) + for (auto j = 0u; j < 3; ++j) + if (fabs(rotation[i][j] - (i == j)) > 1e-12) + SU2_MPI::Error("The first streamwise periodic pair must be purely translational.", CURRENT_FUNCTION); + if (Energy_Equation && Streamwise_Periodic_Temperature && (nMarker_HeatTransfer || nMarker_CHTInterface)) + SU2_MPI::Error("STREAMWISE_PERIODIC_TEMPERATURE requires prescribed MARKER_HEATFLUX walls; " + "MARKER_HEATTRANSFER and MARKER_CHT_INTERFACE are unsupported.", CURRENT_FUNCTION); if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 0) SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION); if (Ref_Inc_NonDim != DIMENSIONAL) diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 0042b62b0f2..66570790729 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -538,7 +538,8 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con /*--- Dot product ---*/ const su2double dot_product = GeometryToolbox::DotProduct(nDim, config->GetPeriodic_Translation(0), Normal); - LinSysRes(iPoint, nDim+1) += scalar_factor*dot_product; + /*--- Normal points into the fluid: q_periodic = q + k grad(T_recovered).n. ---*/ + LinSysRes(iPoint, nDim+1) -= scalar_factor*dot_product; } // if streamwise_periodic break; diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index b6897ef98e1..40921c04af1 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -249,6 +249,12 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, ICOUNT = nPrimVarGrad; JCOUNT = nDim; break; + case PERIODIC_AUXVAR_LS: + COUNT_PER_POINT = nDim*nDim + base_nodes->GetnAuxVar()*nDim; + MPI_TYPE = COMM_TYPE::DOUBLE; + ICOUNT = base_nodes->GetnAuxVar(); + JCOUNT = nDim; + break; case PERIODIC_SOL_LS: case PERIODIC_SOL_ULS: case PERIODIC_SOL_LS_R: @@ -302,6 +308,7 @@ namespace PeriodicCommHelpers { case PERIODIC_PRIM_ULS: return nodes->GetGradient_Primitive(); break; + case PERIODIC_AUXVAR_LS: case PERIODIC_AUXVAR_GG: return nodes->GetAuxVarGradient(); case PERIODIC_SOL_GG: @@ -317,6 +324,7 @@ namespace PeriodicCommHelpers { const su2activematrix& selectField(CVariable* nodes, unsigned short commType) { switch(commType) { + case PERIODIC_AUXVAR_LS: case PERIODIC_AUXVAR_GG: return nodes->GetAuxVar(); case PERIODIC_PRIM_GG: @@ -366,6 +374,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, bool boundary_i, boundary_j; bool weighted = true; + const bool auxiliary = commType == PERIODIC_AUXVAR_LS; + const auto nRotVars = max(max(nVar, nPrimVar), base_nodes->GetnAuxVar()); unsigned short iVar, jVar, iDim; unsigned short nNeighbor = 0; @@ -382,8 +392,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, auto *Und_Lapl = new su2double[nVar]; auto *Sol_Min = new su2double[std::max(nVar, nPrimVarGrad)]; auto *Sol_Max = new su2double[std::max(nVar, nPrimVarGrad)]; - auto *rotPrim_i = new su2double[std::max(nVar, nPrimVar)]; - auto *rotPrim_j = new su2double[std::max(nVar, nPrimVar)]; + auto *rotPrim_i = new su2double[nRotVars]; + auto *rotPrim_j = new su2double[nRotVars]; su2double Sensor_i = 0.0, Sensor_j = 0.0, Pressure_i, Pressure_j; const su2double *Coord_i, *Coord_j; @@ -779,6 +789,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, break; + case PERIODIC_AUXVAR_LS: case PERIODIC_SOL_LS: case PERIODIC_SOL_ULS: case PERIODIC_SOL_LS_R: case PERIODIC_SOL_ULS_R: case PERIODIC_PRIM_LS: case PERIODIC_PRIM_ULS: @@ -821,7 +832,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, for (iVar = 0; iVar < ICOUNT; iVar++) rotPrim_i[iVar] = field(iPoint, iVar); - if (rotate_periodic) { + if (rotate_periodic && !auxiliary) { Rotate(zeros, &field(iPoint,1), &rotPrim_i[1]); } @@ -856,7 +867,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, for (iVar = 0; iVar < ICOUNT; iVar++) rotPrim_j[iVar] = field(jPoint,iVar); - if (rotate_periodic) { + if (rotate_periodic && !auxiliary) { Rotate(zeros, &field(jPoint,1), &rotPrim_j[1]); } @@ -1255,6 +1266,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, break; + case PERIODIC_AUXVAR_LS: case PERIODIC_SOL_LS: case PERIODIC_SOL_ULS: case PERIODIC_SOL_LS_R: case PERIODIC_SOL_ULS_R: case PERIODIC_PRIM_LS: case PERIODIC_PRIM_ULS: @@ -2223,12 +2235,12 @@ void CSolver::SetAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config) void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) { SU2_ZONE_SCOPED - bool weighted = true; + const bool weighted = true; const auto& solution = base_nodes->GetAuxVar(); auto& gradient = base_nodes->GetAuxVarGradient(); auto& rmatrix = base_nodes->GetRmatrix(); - computeGradientsLeastSquares(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, + computeGradientsLeastSquares(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_AUXVAR_LS, *geometry, *config, weighted, solution, 0, base_nodes->GetnAuxVar(), -1, gradient, rmatrix); } diff --git a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp index a981655b39a..53dedec835c 100644 --- a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp @@ -56,6 +56,7 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf, source_cons_jac.resize(nPoint, flamelet_config_options.n_user_scalars) = su2double(0.0); if (flamelet_config_options.preferential_diffusion) { + nAuxVar = FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS) = su2double(0.0); Grad_AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS, nDim, 0.0); } diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp index 11191878554..13e45452679 100644 --- a/UnitTests/Common/CConfig_tests.cpp +++ b/UnitTests/Common/CConfig_tests.cpp @@ -87,3 +87,48 @@ TEST_CASE("INIT_OPTION_INC defaults", "[Config]") { INIT_OPTION_INC::OPERATING_PRESSURE); CHECK(GetInitOptionInc(ideal_gas_options + "INIT_OPTION_INC= DENSITY_INIT\n") == INIT_OPTION_INC::DENSITY_INIT); } + +TEST_CASE("Streamwise periodic Rotation support", "[.StreamwisePeriodicSupportRotation]") { + auto base = base_options; + const std::string custom = "MARKER_CUSTOM= (x_minus, x_plus, z_plus, z_minus)"; + base.replace(base.find(custom), custom.size(), "MARKER_CUSTOM= (z_minus)\nMARKER_SYM= (z_plus)"); + std::stringstream options(base + + "INC_NONDIM= DIMENSIONAL\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP\n" + "INC_ENERGY_EQUATION= YES\nSTREAMWISE_PERIODIC_TEMPERATURE= YES\n" + "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 90,0,0, 1,0,0)\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); +} + +TEST_CASE("Streamwise periodic Zero support", "[.StreamwisePeriodicSupportZero]") { + auto base = base_options; + const std::string custom = "MARKER_CUSTOM= (x_minus, x_plus, z_plus, z_minus)"; + base.replace(base.find(custom), custom.size(), "MARKER_CUSTOM= (z_minus)\nMARKER_SYM= (z_plus)"); + std::stringstream options(base + + "INC_NONDIM= DIMENSIONAL\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP\n" + "INC_ENERGY_EQUATION= YES\nSTREAMWISE_PERIODIC_TEMPERATURE= YES\n" + "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 0,0,0)\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); +} + +TEST_CASE("Streamwise periodic Convection support", "[.StreamwisePeriodicSupportConvection]") { + auto base = base_options; + const std::string custom = "MARKER_CUSTOM= (x_minus, x_plus, z_plus, z_minus)"; + base.replace(base.find(custom), custom.size(), "MARKER_CUSTOM= (z_minus)"); + std::stringstream options( + base + + "INC_NONDIM= DIMENSIONAL\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP\n" + "INC_ENERGY_EQUATION= YES\nSTREAMWISE_PERIODIC_TEMPERATURE= YES\n" + "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)\nMARKER_HEATTRANSFER= (z_plus,1,300)\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); +} + +TEST_CASE("Streamwise periodic CHT support", "[.StreamwisePeriodicSupportCHT]") { + auto base = base_options; + const std::string custom = "MARKER_CUSTOM= (x_minus, x_plus, z_plus, z_minus)"; + base.replace(base.find(custom), custom.size(), "MARKER_CUSTOM= (z_minus)"); + std::stringstream options(base + + "INC_NONDIM= DIMENSIONAL\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP\n" + "INC_ENERGY_EQUATION= YES\nSTREAMWISE_PERIODIC_TEMPERATURE= YES\n" + "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)\nMARKER_CHT_INTERFACE= (z_plus)\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); +} diff --git a/UnitTests/Common/check_periodic_support.py b/UnitTests/Common/check_periodic_support.py new file mode 100644 index 00000000000..5aa6cdf0254 --- /dev/null +++ b/UnitTests/Common/check_periodic_support.py @@ -0,0 +1,50 @@ +#!/usr/bin/env python3 +"""Check periodic configuration errors in isolated processes (MPI errors terminate).""" +import argparse +import subprocess + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("test_driver") + parser.add_argument("--streamwise", action="store_true") + parser.add_argument("--launcher", nargs=argparse.REMAINDER, default=[]) + args = parser.parse_args() + cases = ( + ( + ("Rotation", "must be purely translational"), + ("Zero", "must have a nonzero translation"), + ( + "Convection", + "MARKER_HEATTRANSFER and MARKER_CHT_INTERFACE are unsupported", + ), + ("CHT", "MARKER_HEATTRANSFER and MARKER_CHT_INTERFACE are unsupported"), + ) + if args.streamwise + else ( + ("Continuous", "Continuous adjoints do not implement MARKER_PERIODIC"), + ("Radiation", "RADIATION_MODEL does not implement MARKER_PERIODIC"), + ("Structure", "SOLVER= ELASTICITY does not implement MARKER_PERIODIC"), + ) + ) + prefix = "StreamwisePeriodicSupport" if args.streamwise else "PeriodicSupport" + for tag, message in cases: + result = subprocess.run( + args.launcher + [args.test_driver, "[." + prefix + tag + "]"], + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + text=True, + timeout=30, + check=False, + ) + if result.returncode == 0 or message not in result.stdout: + raise RuntimeError( + tag + + " did not report the expected configuration error:\n" + + result.stdout + ) + print("PASS: " + tag, flush=True) + + +if __name__ == "__main__": + main() diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index 882e363ce5b..3a97294446c 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -26,6 +26,8 @@ */ #include "catch.hpp" +#include "../../SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp" +#include "../../SU2_CFD/include/solvers/CIncNSSolver.hpp" #include "../../Common/include/geometry/CPhysicalGeometry.hpp" #include "../../Common/include/containers/container_decorators.hpp" #include "../../SU2_CFD/include/solvers/CSolver.hpp" @@ -156,16 +158,27 @@ TEST_CASE("LS", "[Gradients]") { testLeastSquares(false); } TEST_CASE("WLS", "[Gradients]") { testLeastSquares(true); } -TEST_CASE("Periodic auxiliary Green-Gauss gradient", "[Gradients][Periodic]") { +TEST_CASE("Periodic auxiliary gradients", "[Gradients][Periodic]") { const bool rotation = GENERATE(false, true); + const bool stretched = GENERATE(false, true); + const auto method = GENERATE(std::string("GREEN_GAUSS"), std::string("WEIGHTED_LEAST_SQUARES")); UnitQuadTestCase field; const string custom = "MARKER_CUSTOM= ( x_minus, x_plus, z_plus, z_minus)"; field.config_options.replace(field.config_options.find(custom), custom.size(), "MARKER_CUSTOM= (z_plus, z_minus)"); field.AddOption(rotation ? "MARKER_PERIODIC= (x_minus,x_plus, 0,0.5,0.5, 90,0,0, 1,0,0)" : "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); field.AddOption("VORTICITY_CONFINEMENT= YES"); + field.AddOption("NUM_METHOD_GRAD= " + method); field.InitConfig(); - field.InitGeometry(); + field.InitGeometry(true); + if (stretched) { + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto x = field.geometry->nodes->GetCoord(i, 0); + field.geometry->nodes->SetCoord(i, 0, x * x); + } + field.geometry->SetControlVolume(field.config.get(), UPDATE); + field.geometry->SetBoundControlVolume(field.config.get(), UPDATE); + } field.geometry->MatchPeriodic(field.config.get(), 1); field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); field.InitSolver(); @@ -177,18 +190,178 @@ TEST_CASE("Periodic auxiliary Green-Gauss gradient", "[Gradients][Periodic]") { const auto* coord = field.geometry->nodes->GetCoord(iPoint); const auto y = coord[1] - 0.5; const auto z = coord[2] - 0.5; - nodes->SetAuxVar(iPoint, 0, rotation ? 2.0 + y * y + z * z : 2.0 + 3.0 * coord[1] - coord[2]); + nodes->SetAuxVar(iPoint, 0, + (rotation ? 2.0 + y * y + z * z : 2.0 + 3.0 * coord[1] - coord[2]) + coord[0] * (1 - coord[0])); + } + if (method == "GREEN_GAUSS") { + SU2_OMP_PARALLEL { solver->SetAuxVar_Gradient_GG(field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + } else { + SU2_OMP_PARALLEL { solver->SetAuxVar_Gradient_LS(field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL } - solver->SetAuxVar_Gradient_GG(field.geometry.get(), field.config.get()); su2double error = 0.0; for (auto iPoint = 0ul; iPoint < field.geometry->GetnPointDomain(); ++iPoint) { const auto* coord = field.geometry->nodes->GetCoord(iPoint); if (rotation && (coord[1] == 0.0 || coord[1] == 1.0 || coord[2] == 0.0 || coord[2] == 1.0)) continue; - const su2double expected[3] = {0.0, rotation ? 2.0 * (coord[1] - 0.5) : 3.0, + const auto index = static_cast( + std::round((stretched ? sqrt(SU2_TYPE::GetValue(coord[0])) : SU2_TYPE::GetValue(coord[0])) * 4)); + auto xcoord = [&](unsigned short k) { + const auto x = 0.25 * k; + return stretched ? x * x : x; + }; + const auto left = xcoord(index == 0 ? 3 : index - 1), right = xcoord(index == 4 ? 1 : index + 1); + const auto hm = coord[0] - (index == 0 ? left - 1 : left); + const auto hp = (index == 4 ? right + 1 : right) - coord[0]; + const auto fm = left * (1 - left), fp = right * (1 - right), fi = coord[0] * (1 - coord[0]); + const auto gx = method == "GREEN_GAUSS" ? (fp - fm) / (hp + hm) : ((fp - fi) / hp + (fi - fm) / hm) / 2; + const su2double expected[3] = {gx, rotation ? 2.0 * (coord[1] - 0.5) : 3.0, rotation ? 2.0 * (coord[2] - 0.5) : -1.0}; for (auto iDim = 0u; iDim < 3; ++iDim) error = max(error, abs(nodes->GetAuxVarGradient(iPoint, 0, iDim) - expected[iDim])); } CHECK(error < 1e-9); } + +TEST_CASE("Flamelet preferential diffusion auxiliary gradients", "[Gradients][Flamelet][Periodic]") { + const bool periodic = GENERATE(false, true); + const auto method = GENERATE(std::string("GREEN_GAUSS"), std::string("WEIGHTED_LEAST_SQUARES")); + UnitQuadTestCase field; + field.SetOption("SOLVER= INC_NAVIER_STOKES"); + field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION"); + field.AddOption("KIND_SCALAR_MODEL= FLAMELET"); + field.AddOption("FLUID_MODEL= FLUID_FLAMELET"); + field.AddOption("INC_DENSITY_MODEL= VARIABLE"); + field.SetOption("VISCOSITY_MODEL= FLAMELET"); + field.AddOption("CONDUCTIVITY_MODEL= FLAMELET"); + field.AddOption("DIFFUSIVITY_MODEL= FLAMELET"); + field.AddOption( + "SPECIES_INIT= (0,300000,0.5)\nCONTROLLING_VARIABLE_NAMES= (ProgressVariable,EnthalpyTot,MixtureFraction)"); + field.AddOption("CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV,NULL,NULL)"); + field.AddOption("PREFERENTIAL_DIFFUSION= YES\nNUM_METHOD_GRAD= " + method); + if (periodic) { + const string custom = "MARKER_CUSTOM= ( x_minus, x_plus, z_plus, z_minus)"; + field.config_options.replace(field.config_options.find(custom), custom.size(), "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); + if (periodic) { + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + } + const su2double initial[] = {0, 300000, 0.5}; + CSpeciesFlameletVariable nodes(initial, field.geometry->GetnPoint(), 3, 3, field.config.get()); + CHECK(nodes.GetnAuxVar() == FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); + /*--- Exercise the shared auxiliary communication with four beta fields and + * three transported scalars, without requiring an unrelated flamelet table. ---*/ + struct AuxiliarySolver : CSolver { + CVariable* nodes; + CVariable* GetBaseClassPointerToNodes() override { return nodes; } + AuxiliarySolver(CVariable& values, CGeometry& geometry) : nodes(&values) { + nDim = 3; + nVar = 3; + nPoint = geometry.GetnPoint(); + nPointDomain = geometry.GetnPointDomain(); + SetBaseClassPointerToNodes(); + } + } solver(nodes, *field.geometry); + if (periodic) { + /*--- The flow solver normally supplies the complete periodic volume. ---*/ + SU2_OMP_PARALLEL { + solver.InitiatePeriodicComms(field.geometry.get(), field.config.get(), 1, PERIODIC_VOLUME); + solver.CompletePeriodicComms(field.geometry.get(), field.config.get(), 1, PERIODIC_VOLUME); + } + END_SU2_OMP_PARALLEL + } + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + for (auto v = 0u; v < FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; ++v) + nodes.SetAuxVar(i, v, (v + 1) * (2 + (periodic ? x[0] * (1 - x[0]) : x[0]) + 3 * x[1] - x[2])); + } + if (method == "GREEN_GAUSS") { + SU2_OMP_PARALLEL { solver.SetAuxVar_Gradient_GG(field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + } else { + SU2_OMP_PARALLEL { solver.SetAuxVar_Gradient_LS(field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + } + su2double error = 0; + for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + const su2double expected[] = {periodic ? (x[0] == 0 || x[0] == 1 ? 0.0 : 1 - 2 * x[0]) : 1.0, 3, -1}; + for (auto v = 0u; v < FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; ++v) + for (auto d = 0u; d < 3; ++d) + error = std::max(error, fabs(nodes.GetAuxVarGradient(i, v, d) - (v + 1) * expected[d])); + } + CHECK(error < 1e-12); +} + +TEST_CASE("Recovered periodic wall heat flux", "[Periodic][Heat]") { + UnitQuadTestCase field; + field.SetOption("SOLVER= INC_NAVIER_STOKES"); + field.AddOption("FLUID_MODEL= CONSTANT_DENSITY"); + field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION"); + field.AddOption("INC_NONDIM= DIMENSIONAL"); + field.AddOption("INC_ENERGY_EQUATION= YES\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP"); + field.AddOption("STREAMWISE_PERIODIC_TEMPERATURE= YES\nINC_VELOCITY_INIT= (1,0,0)"); + const string walls = "MARKER_HEATFLUX= (y_minus, 0.0, y_plus, 0.0)"; + field.config_options.replace(field.config_options.find(walls), walls.size(), "MARKER_HEATFLUX= (y_minus,1,y_plus,3)"); + const string custom = "MARKER_CUSTOM= ( x_minus, x_plus, z_plus, z_minus)"; + field.config_options.replace(field.config_options.find(custom), custom.size(), "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); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + field.geometry->nodes->SetCoord(i, 1, x[1] * (1 + 0.1 * sin(2 * PI_NUMBER * x[0]))); + } + field.geometry->SetControlVolume(field.config.get(), UPDATE); + field.geometry->SetBoundControlVolume(field.config.get(), UPDATE); + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + field.InitSolver(); + auto* solver = static_cast(field.solver[FLOW_SOL]); + field.config->SetGlobalParam(field.config->GetKind_Solver(), RUNTIME_FLOW_SYS); + { + SU2_OMP_PARALLEL { + solver->Preprocessing(field.geometry.get(), field.solver, field.config.get(), MESH_0, 0, RUNTIME_FLOW_SYS, false); + } + END_SU2_OMP_PARALLEL + } + const auto values = solver->GetStreamwisePeriodicValues(); + REQUIRE(fabs(values.Streamwise_Periodic_MassFlow) > 1e-12); + REQUIRE(fabs(values.Streamwise_Periodic_IntegratedHeatFlow) > 1e-12); + auto* nodes = solver->GetNodes(); + su2double error = 0; + unsigned long checked = 0, global = 0; + for (auto m = 0u; m < field.geometry->GetnMarker(); ++m) { + if (field.config->GetMarker_All_TagBound(m) != "y_plus") continue; + solver->LinSysRes = su2double(0); + { + SU2_OMP_PARALLEL { + solver->BC_HeatFlux_Wall(field.geometry.get(), field.solver, nullptr, nullptr, field.config.get(), m); + } + END_SU2_OMP_PARALLEL + } + for (auto v = 0ul; v < field.geometry->GetnVertex(m); ++v) { + const auto i = field.geometry->vertex[m][v]->GetNode(); + if (!field.geometry->nodes->GetDomain(i)) continue; + const auto* normal = field.geometry->vertex[m][v]->GetNormal(); + if (fabs(normal[0]) < 1e-12) continue; + const auto area = GeometryToolbox::Norm(3, normal); + /*--- T = T_periodic + g.x. With an inward area normal, the periodic + * conductive heat flux entering the fluid is q + k g.n. ---*/ + const auto gx = values.Streamwise_Periodic_IntegratedHeatFlow / + (values.Streamwise_Periodic_MassFlow * nodes->GetSpecificHeatCp(i)); + const auto expected = -field.config->GetWall_HeatFlux("y_plus") / field.config->GetHeat_Flux_Ref() * area - + nodes->GetThermalConductivity(i) * gx * normal[0]; + error = std::max(error, fabs(solver->LinSysRes(i, 4) - expected)); + ++checked; + } + } + SU2_MPI::Allreduce(&checked, &global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(global > 0); + CHECK(error < 1e-12); +} diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp index d8da75dc096..d943c1696ee 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()); diff --git a/UnitTests/meson.build b/UnitTests/meson.build index ba3c63afce9..7fa19eb8b3f 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -50,6 +50,9 @@ if get_option('enable-tests') cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args] ) test('Catch2 test driver', test_driver) + test('Streamwise periodic support', python, + args: files('Common/check_periodic_support.py') + [test_driver.full_path(), '--streamwise'], + timeout: 150) test('NEMO convergence field migration', python, args: files('SU2_CFD/output/check_nemo_convergence_migration.py') + [test_driver.full_path()], timeout: 300) From ece527982e03ad32459fe4b7a4b00fbde625eaa0 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 15:53:45 +0200 Subject: [PATCH 09/10] Normalize assignment spacing reported by CodeFactor --- TestCases/parallel_regression.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 6aded1ee082..dfff023d5b4 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -773,7 +773,7 @@ def main(): test_list.append(sp_pipeSlice_3d_dp_hf_tp) # Same with the energy equation, periodic temperature and an integrated wall heat - sp_pipeSlice_3d_dp_ihf_tp = TestCase('sp_pipeSlice_3d_dp_ihf_tp') + sp_pipeSlice_3d_dp_ihf_tp = TestCase('sp_pipeSlice_3d_dp_ihf_tp') sp_pipeSlice_3d_dp_ihf_tp.cfg_dir = "incomp_navierstokes/streamwise_periodic/pipeSlice_3d" sp_pipeSlice_3d_dp_ihf_tp.cfg_file = "sp_pipeSlice_3d_dp_ihf_tp.cfg" sp_pipeSlice_3d_dp_ihf_tp.test_iter = 10 From 95c84b7adbcb3b2cb134bb19fcb85c7078d5f365 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 7 Oct 2026 19:00:46 +0200 Subject: [PATCH 10/10] Promote periodic auxiliary matrix count before multiplication --- 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 40921c04af1..c2f6fb7ddd0 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -250,7 +250,7 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, JCOUNT = nDim; break; case PERIODIC_AUXVAR_LS: - COUNT_PER_POINT = nDim*nDim + base_nodes->GetnAuxVar()*nDim; + COUNT_PER_POINT = static_cast(nDim)*nDim + base_nodes->GetnAuxVar()*nDim; MPI_TYPE = COMM_TYPE::DOUBLE; ICOUNT = base_nodes->GetnAuxVar(); JCOUNT = nDim;