From 583c32afd6953e52afb4cddd83551a6a576a362e Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:09 +0200 Subject: [PATCH 01/12] Add the incompressible grid-motion Jacobian once, with the variable-density terms FDS added the grid-motion term inside a loop over the variables (nVar times). The term is now added once, as the exact derivative with respect to (p, u, h), including d(rho)/dh; the same terms are added to the incompressible JST and Lax-Friedrichs schemes. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- .../src/numerics/flow/convection/centered.cpp | 48 ++++++++++++------- SU2_CFD/src/numerics/flow/convection/fds.cpp | 24 ++++++---- 2 files changed, 46 insertions(+), 26 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/convection/centered.cpp b/SU2_CFD/src/numerics/flow/convection/centered.cpp index b0130279f76f..54fdbdfdd8f9 100644 --- a/SU2_CFD/src/numerics/flow/convection/centered.cpp +++ b/SU2_CFD/src/numerics/flow/convection/centered.cpp @@ -161,17 +161,25 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; /*--- Residual contributions ---*/ - for (iVar = 0; iVar < nVar; iVar++) { + for (iVar = 0; iVar < nVar; iVar++) ProjFlux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); - /*--- Jacobian contributions ---*/ - /*--- Implicit terms ---*/ - if (implicit) { - for (iDim = 0; iDim < nDim+1; iDim++){ - Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added + once. The density depends on h (variable density models) and not on the pressure. ---*/ + if (implicit) { + auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double T, + su2double cp) { + const su2double dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0); + const su2double factor = 0.5*ProjVelocity; + Jac[0][nDim+1] -= factor*dRhodh; + for (iDim = 0; iDim < nDim; iDim++) { + Jac[iDim+1][iDim+1] -= factor*rho; + Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh; } - } + Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh); + }; + addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, Temperature_i, Cp_i); + addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, Temperature_j, Cp_j); } } @@ -401,17 +409,25 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; /*--- Residual contributions ---*/ - for (iVar = 0; iVar < nVar; iVar++) { + for (iVar = 0; iVar < nVar; iVar++) ProjFlux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); - /*--- Jacobian contributions ---*/ - /*--- Implicit terms ---*/ - if (implicit) { - for (iDim = 0; iDim < nDim+1; iDim++){ - Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added + once. The density depends on h (variable density models) and not on the pressure. ---*/ + if (implicit) { + auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double T, + su2double cp) { + const su2double dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0); + const su2double factor = 0.5*ProjVelocity; + Jac[0][nDim+1] -= factor*dRhodh; + for (iDim = 0; iDim < nDim; iDim++) { + Jac[iDim+1][iDim+1] -= factor*rho; + Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh; } - } + Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh); + }; + addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, Temperature_i, Cp_i); + addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, Temperature_j, Cp_j); } } diff --git a/SU2_CFD/src/numerics/flow/convection/fds.cpp b/SU2_CFD/src/numerics/flow/convection/fds.cpp index e806127da203..57c7ed9f6f5a 100644 --- a/SU2_CFD/src/numerics/flow/convection/fds.cpp +++ b/SU2_CFD/src/numerics/flow/convection/fds.cpp @@ -249,19 +249,23 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; /*--- Residual contributions ---*/ - for (iVar = 0; iVar < nVar; iVar++) { + for (iVar = 0; iVar < nVar; iVar++) Flux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); - /*--- Jacobian contributions ---*/ - /*--- Implicit terms ---*/ - if (implicit) { - for (iDim = 0; iDim < nDim; iDim++){ - Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added + once. The density depends on h (variable density models) and not on the pressure. ---*/ + if (implicit) { + auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double dRhodh) { + const su2double factor = 0.5*ProjVelocity; + Jac[0][nDim+1] -= factor*dRhodh; + for (iDim = 0; iDim < nDim; iDim++) { + Jac[iDim+1][iDim+1] -= factor*rho; + Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh; } - Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j; - } + Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh); + }; + addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, dRhodh_i); + addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, dRhodh_j); } } From ca7aee08ccbe731300bbd1e036af55141f00b43d Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:09 +0200 Subject: [PATCH 02/12] Reduce the minimum pressure with MPI_MIN SetRangePressure reduced the local minima with MPI_MAX, so with more than one rank MinP was the largest of the local minima (used by the mass-flow outlet). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 512a99a97dc6..8c6ee1719a01 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2235,7 +2235,7 @@ void CIncEulerSolver::SetRangePressure(CGeometry *geometry, CSolver **solver_con BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { minP = MinP; - SU2_MPI::Allreduce(&minP, &MinP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&minP, &MinP, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm()); maxP = MaxP; SU2_MPI::Allreduce(&maxP, &MaxP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); From b29368df17aee204414315e1d3160ac29bad6c82 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:09 +0200 Subject: [PATCH 03/12] Reject combinations of source terms that share one numerics object The flow solvers build one numerics object per source slot, so two active sources of the same slot overwrote each other (e.g. body force + gravity gave twice the body force). First slot: body force, rotating frame, axisymmetric, gravity, Boussinesq, streamwise periodic. Second slot: radiation, periodic outlet, vorticity confinement. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- Common/src/CConfig.cpp | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 034352b92aa4..c5fef6547554 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4598,6 +4598,36 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i Rotating_Frame = (Kind_GridMovement == ROTATING_FRAME); + /*--- The finite volume flow solvers apply these source terms with a single numerics object, which is built + for only one of them: combinations would apply one source several times and ignore the others. Only these + solvers (and their discrete adjoints) are checked: e.g. elasticity uses BODY_FORCE and GRAVITY_FORCE together. ---*/ + + if (GetFluidProblem() && !GetNEMOProblem()) { + const bool inc = (Kind_Regime == ENUM_REGIME::INCOMPRESSIBLE); + const int nSources = int(Body_Force) + int(Rotating_Frame) + int(Axisymmetric) + int(GravityForce && !inc) + + int(inc && Kind_DensityModel == INC_DENSITYMODEL::BOUSSINESQ) + + int(inc && Kind_Streamwise_Periodic != ENUM_STREAMWISE_PERIODIC::NONE); + if (nSources > 1) { + SU2_MPI::Error("Only one of BODY_FORCE, GRID_MOVEMENT= ROTATING_FRAME, AXISYMMETRIC, GRAVITY_FORCE (compressible), " + "INC_DENSITY_MODEL= BOUSSINESQ and KIND_STREAMWISE_PERIODIC (incompressible) can be used at a time.", + CURRENT_FUNCTION); + } + } + + /*--- The second source object, for radiation, the energy correction of incompressible streamwise periodic flow + (without periodic temperature) or vorticity confinement, is also built for only one of them. ---*/ + + if (GetFluidProblem() && !GetNEMOProblem()) { + const bool inc = (Kind_Regime == ENUM_REGIME::INCOMPRESSIBLE); + const bool periodicOutlet = inc && Kind_Streamwise_Periodic != ENUM_STREAMWISE_PERIODIC::NONE && Energy_Equation && + !Streamwise_Periodic_Temperature; + if (int(Radiation) + int(periodicOutlet) + int(VorticityConfinement) > 1) { + SU2_MPI::Error("Only one of RADIATION, the energy correction of KIND_STREAMWISE_PERIODIC (incompressible, " + "STREAMWISE_PERIODIC_TEMPERATURE= NO) and VORTICITY_CONFINEMENT can be used at a time.", + CURRENT_FUNCTION); + } + } + /*--- In case the grid movement parameters have not been declared in the config file, set them equal to zero for safety. Also check to make sure that for each option, a value has been declared for each moving marker. ---*/ From df8f9787cdea64283d369dfd773793b10c63557c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:09 +0200 Subject: [PATCH 04/12] Read the Riemann boundary kind in BC_TurboRiemann It read the Giles kind of the marker (copied from BC_Giles), which reads past the end of the Giles list and crashes when the case has no Giles marker. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/src/solvers/CEulerSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 8259e1630520..15658e5a7459 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -5800,7 +5800,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim]; su2double donorAverages[5] = {0.0}; - switch (config->GetKind_Data_Giles(Marker_Tag)){ + switch (config->GetKind_Data_Riemann(Marker_Tag)){ case MIXING_IN: case MIXING_IN_1D: case MIXING_OUT: case MIXING_OUT_1D: for (auto mixVar = 0u; mixVar < 5; mixVar++) donorAverages[mixVar] = GetMixingState(val_marker, iSpan, mixVar); break; From f3bcaac94b17d55ed5217ee46f350089b794a26b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:09 +0200 Subject: [PATCH 05/12] Check the number of variables in the Green-Gauss gradient The boundary flux array had 20 entries and no bound check (NEMO with 10 or more species in 3D). Use 32, as the limiters, and stop with an error above it. The error of both checks is raised by the master thread only, since these functions run inside OpenMP parallel regions (MPI_THREAD_FUNNELED). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp | 9 ++++++++- SU2_CFD/include/limiters/computeLimiters_impl.hpp | 6 +++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp index f60c0422f0cc..89eb028c6aeb 100644 --- a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp +++ b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp @@ -69,7 +69,14 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER const auto chunkSize = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_CHUNK); #endif - static constexpr size_t MAXNVAR = 20; + static constexpr size_t MAXNVAR = 32; + + /*--- Called by all threads of a parallel region: only the master may call MPI (MPI_THREAD_FUNNELED). ---*/ + if (varEnd > MAXNVAR) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + SU2_MPI::Error("Number of variables is too large, increase MAXNVAR.", CURRENT_FUNCTION); + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } /*--- For each (non-halo) volume integrate over its faces (edges). ---*/ diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp index 7d8beafb03e0..0cdc60207a27 100644 --- a/SU2_CFD/include/limiters/computeLimiters_impl.hpp +++ b/SU2_CFD/include/limiters/computeLimiters_impl.hpp @@ -80,8 +80,12 @@ void computeLimiters_impl(CSolver* solver, { constexpr size_t MAXNVAR = 32; - if (varEnd > MAXNVAR) + /*--- Called by all threads of a parallel region: only the master may call MPI (MPI_THREAD_FUNNELED). ---*/ + if (varEnd > MAXNVAR) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS SU2_MPI::Error("Number of variables is too large, increase MAXNVAR.", CURRENT_FUNCTION); + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } const size_t nPointDomain = geometry.GetnPointDomain(); const size_t nPoint = geometry.GetnPoint(); From 99d861a2e510e9ea8edf735f7d0a649431b50a03 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:10 +0200 Subject: [PATCH 06/12] Clamp MGLEVEL to 10 The multigrid arrays hold at most 10 coarse levels (MAX_MG_LEVELS); MGLEVEL= 12 with MG_MIN_MESHSIZE= 1 crashed with stack smashing. Clamp with a warning, MG_MIN_MESHSIZE can still reduce the levels. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- Common/src/CConfig.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c5fef6547554..3a5bfee8f713 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4987,6 +4987,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (Restart || ((Kind_MGCycle == MG_CYCLE::FULL) && ContinuousAdjoint)) Kind_MGCycle = MG_CYCLE::V; + /*--- The multigrid arrays hold at most 10 coarse levels (CMultiGridIntegration::MAX_MG_LEVELS). Clamp instead of + stopping, MG_MIN_MESHSIZE can still reduce the number of levels during the agglomeration. ---*/ + if (nMGLevels > 10) { + if (rank == MASTER_NODE) cout << "WARNING: MGLEVEL > 10 is not supported, using MGLEVEL= 10." << endl; + nMGLevels = 10; + } + FinestMesh = MESH_0; if (Kind_MGCycle == MG_CYCLE::FULL) FinestMesh = nMGLevels; From 536c8d8578fe6128b0492963357474f80a2dd62a Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:09:20 +0200 Subject: [PATCH 07/12] Use 0.5 instead of the integer division 1/2 in the axisymmetric source Jacobian The term 1/2*(Gamma-1) of d(v (rho E + p))/d(rho v) was 0 (implicit Jacobian only). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index e3797c7c5a5f..87784e73f39f 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -138,7 +138,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi jacobian[3][0] = -Gamma*U_i[2]*U_i[3]/(U_i[0]*U_i[0]) + (Gamma-1)*U_i[2]*(U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]*U_i[0]); jacobian[3][1] = -(Gamma-1)*U_i[2]*U_i[1]/(U_i[0]*U_i[0]); - jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 1/2*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) ); + jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 0.5*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) ); jacobian[3][3] = Gamma*U_i[2]/U_i[0]; for (iVar=0; iVar < nVar; iVar++) From c333a18ebd6c26c584e24088fedfe9b1c9364af5 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 29 Sep 2026 14:10:40 +0200 Subject: [PATCH 08/12] Check the MUSCL face states before evaluating the fluid model The consistent extrapolation and the low-Mach correction called the fluid model with the reconstructed pressure and density before checking that they are positive (possible NaN or fluid-model errors with real-gas models at start-up). Check them first, also rejecting zero and NaN, and check again after the low-Mach correction. Valid states give the same result as before. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/src/solvers/CEulerSolver.cpp | 42 ++++++++++++++++++---------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 15658e5a7459..6b8784f228aa 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1947,27 +1947,39 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain } - /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. ---*/ + /*--- Check for non-physical solutions after reconstruction. If found, use the + cell-average value of the solution. This is a locally 1st order approximation, + which is typically only active during the start-up of a calculation. + Pressure and density are checked before the fluid model is evaluated with them. ---*/ - if (!ideal_gas || low_mach_corr) { - ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i); - ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j); - } + auto nonPositive = [](su2double value) { return !(value > 0.0) || !std::isfinite(SU2_TYPE::GetValue(value)); }; + bool neg_pres_or_rho_i = nonPositive(Primitive_i[prim_idx.Pressure()]) || nonPositive(Primitive_i[prim_idx.Density()]); + bool neg_pres_or_rho_j = nonPositive(Primitive_j[prim_idx.Pressure()]) || nonPositive(Primitive_j[prim_idx.Density()]); - /*--- Low-Mach number correction. ---*/ + if (!neg_pres_or_rho_i && !neg_pres_or_rho_j) { - if (low_mach_corr) { - LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j); - } + /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. ---*/ - /*--- Check for non-physical solutions after reconstruction. If found, use the - cell-average value of the solution. This is a locally 1st order approximation, - which is typically only active during the start-up of a calculation. ---*/ + if (!ideal_gas || low_mach_corr) { + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i); + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j); + } + + /*--- Low-Mach number correction. ---*/ + + if (low_mach_corr) { + LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j); + } + + /*--- The correction changes the face states, check them again. ---*/ - bool neg_pres_or_rho_i = (Primitive_i[prim_idx.Pressure()] < 0.0) || (Primitive_i[prim_idx.Density()] < 0.0); - bool neg_pres_or_rho_j = (Primitive_j[prim_idx.Pressure()] < 0.0) || (Primitive_j[prim_idx.Density()] < 0.0); + neg_pres_or_rho_i = nonPositive(Primitive_i[prim_idx.Pressure()]) || nonPositive(Primitive_i[prim_idx.Density()]); + neg_pres_or_rho_j = nonPositive(Primitive_j[prim_idx.Pressure()]) || nonPositive(Primitive_j[prim_idx.Density()]); + } - su2double R = sqrt(fabs(Primitive_j[prim_idx.Density()]/Primitive_i[prim_idx.Density()])); + su2double R = 1.0; + if (!neg_pres_or_rho_i && !neg_pres_or_rho_j) + R = sqrt(Primitive_j[prim_idx.Density()]/Primitive_i[prim_idx.Density()]); su2double sq_vel = 0.0; for (auto iDim = 0u; iDim < nDim; iDim++) { su2double RoeVelocity = (R * Primitive_j[iDim + prim_idx.Velocity()] + From e0d3ee8d72598767227e62a9be5834a4554532b6 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 21:19:39 +0200 Subject: [PATCH 09/12] Correct compressible moving-wall projection and multigrid coupling --- .../include/solvers/CFVMFlowSolverBase.inl | 45 ++++- .../SU2_CFD/solvers/moving_wall_tests.cpp | 182 ++++++++++++++++++ UnitTests/meson.build | 1 + 3 files changed, 218 insertions(+), 10 deletions(-) create mode 100644 UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index cb44757c6280..650b26b91e7f 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1155,6 +1155,7 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || (config->GetKind_FluidModel() == IDEAL_GAS); const auto iVel = prim_idx.Velocity(); + const bool compressibleProjection = FlowRegime == ENUM_REGIME::COMPRESSIBLE && iVel == 1 && nVar == nDim + 2; /*--- Blazek chapter 8.: * The components of the momentum residual normal to the symmetry plane are zeroed out. @@ -1184,8 +1185,8 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area; } - /*--- Energy terms due to grid movement (aka work of pressure forces). ---*/ - if (dynamic_grid) { + /*--- The classical compressible projection below already supplies wall pressure work. ---*/ + if (dynamic_grid && !compressibleProjection) { su2double* V_reflected = GetCharacPrimVar(val_marker, iVertex); conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), @@ -1217,8 +1218,7 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve /*--- Compute the residual using an upwind scheme. ---*/ auto residual = conv_numerics->ComputeResidual(config); - /*--- Use just the energy fluxes to update the residual, adding the others would - * increase numerical diffusion which we wish to avoid if possible. ---*/ + /*--- Energy fluxes for the other conservative layouts. ---*/ for (auto iVar = iVel + nDim; iVar < nVar; iVar++) { LinSysRes(iPoint, iVar) += residual.residual[iVar]; } @@ -1243,14 +1243,16 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve su2double* solutionOld = nodes->GetSolution_Old(iPoint); su2double gridVel[MAXNDIM] = {}; + su2double gridNormalVel = 0.0; if (dynamic_grid) { for (auto iDim = 0u; iDim < nDim; iDim++) { gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim]; } + gridNormalVel = GeometryToolbox::DotProduct(nDim, gridVel, UnitNormal); if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) { for(auto iDim = 0u; iDim < nDim; iDim++) { /*--- Multiply by density since we are correcting conservative variables. ---*/ - gridVel[iDim] *= nodes->GetDensity(iPoint); + gridVel[iDim] *= compressibleProjection ? solutionOld[0] : nodes->GetDensity(iPoint); } } } @@ -1261,20 +1263,27 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve for (auto iDim = 0u; iDim < nDim; iDim++) { solutionOld[iVel + iDim] -= vp * UnitNormal[iDim]; } + if (compressibleProjection) { + solutionOld[iVel + nDim] -= gridNormalVel * vp; + } /*--- Keep only the tangential part of the momentum residuals. ---*/ su2double normalRes = 0.0; for (auto iDim = 0u; iDim < nDim; iDim++) { normalRes += LinSysRes(iPoint, iVel + iDim) * UnitNormal[iDim]; } + if (compressibleProjection) { + /*--- Project in the wall frame: normal momentum follows density, with consistent pressure work. ---*/ + normalRes -= gridNormalVel * LinSysRes(iPoint, 0); + LinSysRes(iPoint, iVel + nDim) -= gridNormalVel * normalRes; + } for (auto iDim = 0u; iDim < nDim; iDim++) { LinSysRes(iPoint, iVel + iDim) -= normalRes * UnitNormal[iDim]; } /*--- Jacobian contribution for implicit integration. ---*/ if (implicit) { - /*--- Modify the Jacobians according to the modification of the residual - * J_new = (I - n * n^T) * J where n = {0, nx, ny, nz, 0, ...} ---*/ + /*--- Apply the same wall-frame projection to the Jacobians. ---*/ su2double mat[MAXNVAR * MAXNVAR] = {}; for (auto iVar = 0u; iVar < nVar; iVar++) @@ -1282,6 +1291,13 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve for (auto iDim = 0u; iDim < nDim; iDim++) for (auto jDim = 0u; jDim < nDim; jDim++) mat[(iDim + iVel) * nVar + jDim + iVel] -= UnitNormal[iDim] * UnitNormal[jDim]; + if (compressibleProjection) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + mat[(iDim + iVel) * nVar] += gridNormalVel * UnitNormal[iDim]; + mat[(iVel + nDim) * nVar + iVel + iDim] -= gridNormalVel * UnitNormal[iDim]; + } + mat[(iVel + nDim) * nVar] += gridNormalVel * gridNormalVel; + } auto ModifyJacobian = [&](const unsigned long jPoint) { su2double jac[MAXNVAR * MAXNVAR], newJac[MAXNVAR * MAXNVAR]; @@ -1306,6 +1322,10 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve for (auto iDim = 0u; iDim < nDim; iDim++) { normalRes += Res_TruncError[iVel + iDim] * UnitNormal[iDim]; } + if (compressibleProjection) { + normalRes -= gridNormalVel * Res_TruncError[0]; + Res_TruncError[iVel + nDim] -= gridNormalVel * normalRes; + } for (auto iDim = 0u; iDim < nDim; iDim++) { Res_TruncError[iVel + iDim] -= normalRes * UnitNormal[iDim]; } @@ -2974,6 +2994,7 @@ void CFVMFlowSolverBase::MultigridProjectEulerWall(CGeometry* geo bool use_solution_old) { const auto iVel = prim_idx.Velocity(); const auto nDim = geometry->GetnDim(); + const bool compressibleProjection = FlowRegime == ENUM_REGIME::COMPRESSIBLE && iVel == 1 && nVar == nDim + 2; for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetMarker_All_KindBC(iMarker) != EULER_WALL) continue; @@ -3004,12 +3025,15 @@ void CFVMFlowSolverBase::MultigridProjectEulerWall(CGeometry* geo * For dynamic grids subtract the grid velocity to enforce (v - v_grid).n = 0, * multiplying by density for compressible flow (conservative variables). ---*/ su2double gridVel[MAXNDIM] = {}; - if (dynamic_grid && !use_solution_old) { + su2double gridNormalVel = 0.0; + if (dynamic_grid && (compressibleProjection || !use_solution_old)) { for (auto iDim = 0u; iDim < nDim; iDim++) gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim]; + gridNormalVel = GeometryToolbox::DotProduct(nDim, gridVel, UnitNormal); if constexpr (FlowRegime == ENUM_REGIME::COMPRESSIBLE) { - for (auto iDim = 0u; iDim < nDim; iDim++) - gridVel[iDim] *= nodes->GetDensity(iPoint); + /*--- Corrections carry their own density increment, rather than the current solution density. ---*/ + const auto density = compressibleProjection ? sol[0] : nodes->GetDensity(iPoint); + for (auto iDim = 0u; iDim < nDim; iDim++) gridVel[iDim] *= density; } } @@ -3019,6 +3043,7 @@ void CFVMFlowSolverBase::MultigridProjectEulerWall(CGeometry* geo /*--- Project to tangent plane. ---*/ for (auto iDim = 0u; iDim < nDim; iDim++) sol[iVel + iDim] -= momentum_n * UnitNormal[iDim]; + if (compressibleProjection) sol[iVel + nDim] -= gridNormalVel * momentum_n; } END_SU2_OMP_FOR } diff --git a/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp b/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp new file mode 100644 index 000000000000..1baf0c271428 --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp @@ -0,0 +1,182 @@ +/*! + * \file moving_wall_tests.cpp + * \brief Regression for compressible wall projection in fixed and uniformly translating frames. + * \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 +#include +#include "../../UnitQuadTestCase.hpp" +#include "../../../SU2_CFD/include/numerics/flow/convection/hllc.hpp" +#include "../../../SU2_CFD/include/fluid/CPengRobinson.hpp" + +TEST_CASE("Compressible wall projection commutes with a Galilean boost", "[HLLC-wall]") { + for (const auto nDim : {2u, 3u}) { + for (const bool realGas : {false, true}) { + UnitQuadTestCase fixed, moving; + const std::string markers = + nDim == 2 ? "(x_minus,x_plus,y_minus,y_plus)" : "(x_minus,x_plus,y_minus,y_plus,z_minus,z_plus)"; + const std::string options = "SOLVER= EULER\nMESH_FORMAT= " + std::string(nDim == 2 ? "RECTANGLE" : "BOX") + + "\nMESH_BOX_SIZE= 5,5,5\nREF_DIMENSIONALIZATION= DIMENSIONAL\n" + "MACH_NUMBER= 0.1\nFREESTREAM_OPTION= DENSITY_FS\nFREESTREAM_DENSITY= 230\n" + "FREESTREAM_PRESSURE= 1500112\nGAS_CONSTANT= 35.152\n" + "CONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= EULER_IMPLICIT\n" + "GRID_MOVEMENT= STEADY_TRANSLATION\nTRANSLATION_RATE= 0,0,0\nMARKER_EULER= " + + markers + "\nMARKER_MONITORING= (x_minus)\n" + + (realGas ? "FLUID_MODEL= PR_GAS\nGAMMA_VALUE= 1.0125\nCRITICAL_PRESSURE= " + "1415200\nCRITICAL_TEMPERATURE= 564.1\nACENTRIC_FACTOR= 0.529\n" + : "FLUID_MODEL= IDEAL_GAS\nGAMMA_VALUE= 1.4\n"); + fixed.config_options = moving.config_options = options; + fixed.InitConfig(); + fixed.InitGeometry(); + fixed.InitSolver(); + moving.InitConfig(); + moving.InitGeometry(); + moving.InitSolver(); + const auto nVar = nDim + 2; + const su2double boost[3] = {20, -7, 3}; + const su2double velocity[3] = {11, 4, -2}; + su2double B[5][5] = {}, T[5][5] = {}, baseJac[5][5] = {}, boostedJac[5][5] = {}; + su2double c2 = 0; + for (auto i = 0u; i < nVar; ++i) { + B[i][i] = T[i][i] = 1; + for (auto j = 0u; j < nVar; ++j) baseJac[i][j] = 1 + i + 2 * j; + } + for (auto i = 0u; i < nDim; ++i) { + c2 += boost[i] * boost[i]; + B[i + 1][0] = boost[i]; + T[i + 1][0] = -boost[i]; + B[nVar - 1][i + 1] = boost[i]; + T[nVar - 1][i + 1] = -boost[i]; + } + B[nVar - 1][0] = T[nVar - 1][0] = 0.5 * c2; + for (auto i = 0u; i < nVar; ++i) + for (auto j = 0u; j < nVar; ++j) + for (auto k = 0u; k < nVar; ++k) + for (auto l = 0u; l < nVar; ++l) boostedJac[i][j] += B[i][k] * baseJac[k][l] * T[l][j]; + std::unique_ptr fluid = + realGas ? std::unique_ptr(new CPengRobinson(1.0125, 35.152, 1415200, 564.1, 0.529)) + : std::unique_ptr(new CIdealGas(1.4, 35.152)); + fluid->SetTDState_Prho(1500112, 230); + std::array baseState{}; + baseState[0] = 230; + su2double speed2 = 0; + for (auto i = 0u; i < nDim; ++i) { + baseState[i + 1] = 230 * velocity[i]; + speed2 += velocity[i] * velocity[i]; + } + baseState[nVar - 1] = 230 * (fluid->GetStaticEnergy() + 0.5 * speed2); + const su2double baseResidual[5] = {0.3, -0.5, 0.2, 0.7, -0.1}; + for (auto frame = 0u; frame < 2; ++frame) { + auto& test = frame ? moving : fixed; + auto* solver = test.solver[FLOW_SOL]; + auto* nodes = static_cast(solver->GetNodes()); + for (auto p = 0ul; p < test.geometry->GetnPoint(); ++p) { + su2double state[5] = {}, residual[5] = {}; + for (auto i = 0u; i < nVar; ++i) + for (auto j = 0u; j < nVar; ++j) { + state[i] += (frame ? B[i][j] : (i == j)) * baseState[j]; + residual[i] += (frame ? B[i][j] : (i == j)) * baseResidual[j]; + } + nodes->SetSolution(p, state); + nodes->SetSolution_Old(p, state); + REQUIRE(nodes->SetPrimVar(p, fluid.get())); + nodes->SetSecondaryVar(p, fluid.get()); + for (auto i = 0u; i < nDim; ++i) test.geometry->nodes->SetGridVel(p, i, frame ? boost[i] : 0); + for (auto i = 0u; i < nVar; ++i) { + solver->LinSysRes(p, i) = residual[i]; + nodes->GetResTruncError(p)[i] = residual[i]; + } + su2double block[25] = {}; + for (auto i = 0u; i < nVar; ++i) + for (auto j = 0u; j < nVar; ++j) block[i * nVar + j] = frame ? boostedJac[i][j] : baseJac[i][j]; + solver->Jacobian.SetBlock(p, p, block); + for (auto k = 0u; k < test.geometry->nodes->GetnPoint(p); ++k) + solver->Jacobian.SetBlock(p, test.geometry->nodes->GetPoint(p, k), block); + } + std::unique_ptr numerics = + realGas ? std::unique_ptr(new CUpwGeneralHLLC_Flow(nDim, nVar, test.config.get())) + : std::unique_ptr(new CUpwHLLC_Flow(nDim, nVar, test.config.get())); + for (auto marker = 0u; marker < test.config->GetnMarker_All(); ++marker) + if (test.config->GetMarker_All_TagBound(marker) == "x_minus") + solver->BC_Sym_Plane(test.geometry.get(), test.solver, numerics.get(), nullptr, test.config.get(), marker); + } + for (auto p = 0ul; p < fixed.geometry->GetnPoint(); ++p) { + auto* f = fixed.solver[FLOW_SOL]; + auto* m = moving.solver[FLOW_SOL]; + for (auto i = 0u; i < nVar; ++i) { + su2double state = 0, residual = 0, truncation = 0; + for (auto j = 0u; j < nVar; ++j) { + state += B[i][j] * f->GetNodes()->GetSolution_Old(p)[j]; + residual += B[i][j] * f->LinSysRes(p, j); + truncation += B[i][j] * f->GetNodes()->GetResTruncError(p)[j]; + } + CHECK(m->GetNodes()->GetSolution_Old(p)[i] == Approx(state).margin(1e-6)); + CHECK(m->LinSysRes(p, i) == Approx(residual).margin(1e-8)); + CHECK(m->GetNodes()->GetResTruncError(p)[i] == Approx(truncation).margin(1e-8)); + for (auto block = 0u; block <= fixed.geometry->nodes->GetnPoint(p); ++block) { + const auto neighbor = block ? fixed.geometry->nodes->GetPoint(p, block - 1) : p; + for (auto j = 0u; j < nVar; ++j) { + su2double expected = 0; + for (auto k = 0u; k < nVar; ++k) + for (auto l = 0u; l < nVar; ++l) + expected += B[i][k] * f->Jacobian.GetBlock(p, neighbor)[k * nVar + l] * T[l][j]; + CHECK(m->Jacobian.GetBlock(p, neighbor)[i * nVar + j] == Approx(expected).margin(1e-5)); + } + } + } + } + + for (const bool correction : {false, true}) { + const su2double delta[5] = {0.3, -0.5, 0.2, 0.7, -0.1}; + for (auto frame = 0u; frame < 2; ++frame) { + auto& test = frame ? moving : fixed; + auto* solver = test.solver[FLOW_SOL]; + for (auto p = 0ul; p < test.geometry->GetnPoint(); ++p) { + su2double state[5] = {}; + for (auto i = 0u; i < nVar; ++i) + for (auto j = 0u; j < nVar; ++j) + state[i] += (frame ? B[i][j] : (i == j)) * (correction ? delta[j] : baseState[j]); + if (correction) + solver->GetNodes()->SetSolution_Old(p, state); + else + solver->GetNodes()->SetSolution(p, state); + } + solver->MultigridProjectEulerWall(test.geometry.get(), test.config.get(), correction); + } + for (auto p = 0ul; p < fixed.geometry->GetnPoint(); ++p) { + const auto* f = correction ? fixed.solver[FLOW_SOL]->GetNodes()->GetSolution_Old(p) + : fixed.solver[FLOW_SOL]->GetNodes()->GetSolution(p); + const auto* m = correction ? moving.solver[FLOW_SOL]->GetNodes()->GetSolution_Old(p) + : moving.solver[FLOW_SOL]->GetNodes()->GetSolution(p); + for (auto i = 0u; i < nVar; ++i) { + su2double expected = 0; + for (auto j = 0u; j < nVar; ++j) expected += B[i][j] * f[j]; + CHECK(m[i] == Approx(expected).margin(1e-6)); + } + } + } + } + } +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index bb5448765af6..46a9c9b3cbbc 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -22,6 +22,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'SU2_CFD/output/COutput_convergence_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp', + 'SU2_CFD/solvers/moving_wall_tests.cpp', 'SU2_CFD/nemo_viscous_assembly.cpp', 'SU2_CFD/windowing.cpp', 'Common/toolboxes/random_toolbox_tests.cpp', From a950f3d3cbea0501415a76aa2484486d9278b68b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 22:29:00 +0200 Subject: [PATCH 10/12] Reject unsupported multigrid levels and align centered density handling --- Common/src/CConfig.cpp | 6 ++---- .../src/numerics/flow/convection/centered.cpp | 8 ++++---- UnitTests/Common/CConfig_tests.cpp | 14 ++++++++++++++ UnitTests/Common/check_multigrid_limits.py | 16 ++++++++++++++++ UnitTests/meson.build | 2 ++ 5 files changed, 38 insertions(+), 8 deletions(-) create mode 100644 UnitTests/Common/check_multigrid_limits.py diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index ae6f75fbc969..be485288aaf0 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4994,11 +4994,9 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (Restart || ((Kind_MGCycle == MG_CYCLE::FULL) && ContinuousAdjoint)) Kind_MGCycle = MG_CYCLE::V; - /*--- The multigrid arrays hold at most 10 coarse levels (CMultiGridIntegration::MAX_MG_LEVELS). Clamp instead of - stopping, MG_MIN_MESHSIZE can still reduce the number of levels during the agglomeration. ---*/ + /*--- The multigrid arrays hold at most 10 coarse levels (CMultiGridIntegration::MAX_MG_LEVELS). ---*/ if (nMGLevels > 10) { - if (rank == MASTER_NODE) cout << "WARNING: MGLEVEL > 10 is not supported, using MGLEVEL= 10." << endl; - nMGLevels = 10; + SU2_MPI::Error("MGLEVEL must be between 0 and 10; larger multigrid hierarchies are not supported.", CURRENT_FUNCTION); } FinestMesh = MESH_0; diff --git a/SU2_CFD/src/numerics/flow/convection/centered.cpp b/SU2_CFD/src/numerics/flow/convection/centered.cpp index 54fdbdfdd8f9..713ca2e86662 100644 --- a/SU2_CFD/src/numerics/flow/convection/centered.cpp +++ b/SU2_CFD/src/numerics/flow/convection/centered.cpp @@ -31,7 +31,7 @@ CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - variable_density = ((config->GetKind_DensityModel() == INC_DENSITYMODEL::VARIABLE) || (config->GetKind_DensityModel() == INC_DENSITYMODEL::FLAMELET)); + variable_density = config->GetVariable_Density_Model(); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); fix_factor = config->GetCent_Inc_Jac_Fix_Factor(); @@ -169,7 +169,7 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi if (implicit) { auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double T, su2double cp) { - const su2double dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0); + const auto dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0); const su2double factor = 0.5*ProjVelocity; Jac[0][nDim+1] -= factor*dRhodh; for (iDim = 0; iDim < nDim; iDim++) { @@ -258,7 +258,7 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi CCentJSTInc_Flow::CCentJSTInc_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - variable_density = (config->GetVariable_Density_Model()); + variable_density = config->GetVariable_Density_Model(); energy = config->GetEnergy_Equation(); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); @@ -417,7 +417,7 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi if (implicit) { auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double T, su2double cp) { - const su2double dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0); + const auto dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0); const su2double factor = 0.5*ProjVelocity; Jac[0][nDim+1] -= factor*dRhodh; for (iDim = 0; iDim < nDim; iDim++) { diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp index 111918785545..20a34724fadc 100644 --- a/UnitTests/Common/CConfig_tests.cpp +++ b/UnitTests/Common/CConfig_tests.cpp @@ -87,3 +87,17 @@ 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("Supported multigrid levels are retained", "[Config][MultigridLimits]") { + for (const auto levels : {0, 10}) { + std::stringstream options(base_options + "MGLEVEL= " + std::to_string(levels) + "\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + CHECK(config.GetnMGLevels() == levels); + } +} + +TEST_CASE("Excessive multigrid levels are rejected", "[.MultigridLimitError]") { + std::stringstream options(base_options + "MGLEVEL= 11\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + FAIL("An unsupported multigrid hierarchy was accepted."); +} diff --git a/UnitTests/Common/check_multigrid_limits.py b/UnitTests/Common/check_multigrid_limits.py new file mode 100644 index 000000000000..5070fdaed3cc --- /dev/null +++ b/UnitTests/Common/check_multigrid_limits.py @@ -0,0 +1,16 @@ +#!/usr/bin/env python3 +"""Check the fatal configuration diagnostic in its own process.""" +import subprocess +import sys + +result = subprocess.run( + [sys.argv[1], "[.MultigridLimitError]"], capture_output=True, text=True, timeout=30 +) +expected = "MGLEVEL must be between 0 and 10" +if result.returncode == 0 or expected not in result.stdout + result.stderr: + raise RuntimeError( + "Unsupported MGLEVEL did not fail with its diagnostic:\n" + + result.stdout + + result.stderr + ) +print("PASS: MGLEVEL=11 is rejected with its supported range") diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 46a9c9b3cbbc..de7945480f08 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -53,6 +53,8 @@ if get_option('enable-tests') cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args] ) test('Catch2 test driver', test_driver) + test('Multigrid configuration limits', python, + args: files('Common/check_multigrid_limits.py') + [test_driver.full_path()]) test('NEMO convergence field migration', python, args: files('SU2_CFD/output/check_nemo_convergence_migration.py') + [test_driver.full_path()], timeout: 300) From 27fc17589a63258d16db9589bc58df540b313355 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 22:28:44 +0200 Subject: [PATCH 11/12] Separate the HLLC Jacobian option and keep verification access private --- Common/include/CConfig.hpp | 4 +++ Common/src/CConfig.cpp | 4 ++- .../include/numerics/flow/convection/hllc.hpp | 30 +++---------------- SU2_CFD/src/numerics/flow/convection/hllc.cpp | 12 +++----- UnitTests/Common/CConfig_tests.cpp | 16 ++++++++++ .../SU2_CFD/numerics/hllc_jacobian_tests.cpp | 19 ++++++++---- config_template.cfg | 5 +++- 7 files changed, 49 insertions(+), 41 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 9d4920c871b0..183ecde303af 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -607,6 +607,7 @@ class CConfig { MUSCL_Kappa_Heat, /*!< \brief Blending coefficient for U-MUSCL scheme for the (fvm) heat equation.*/ MUSCL_Kappa_AdjFlow, /*!< \brief Blending coefficient for U-MUSCL scheme for the adj flow equations.*/ MUSCL_Kappa_Species; /*!< \brief Blending coefficient for U-MUSCL scheme for the species equations.*/ + bool Use_Exact_HLLC_Jacobians; /*!< \brief Exact HLLC Jacobians with wave speeds held fixed. */ bool Use_Accurate_Jacobians; /*!< \brief Use numerically computed Jacobians for AUSM+up(2) and SLAU(2). */ bool Use_Accurate_Turb_Jacobians; /*!< \brief Use numerically computed Jacobians for standard SA turbulence model. */ bool EulerPersson; /*!< \brief Boolean to determine whether this is an Euler simulation with Persson shock capturing. */ @@ -4800,6 +4801,9 @@ class CConfig { */ bool GetUse_Accurate_Jacobians(void) const { return Use_Accurate_Jacobians; } + /*!\brief Get whether to use exact HLLC Jacobians with wave speeds held fixed. */ + bool GetUse_Exact_HLLC_Jacobians(void) const { return Use_Exact_HLLC_Jacobians; } + /*! * \brief Get whether to "Use Accurate Jacobians" for Standard SA turbulence model. * \return yes/no. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index be485288aaf0..81fb736ce20d 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2176,8 +2176,10 @@ void CConfig::SetConfig_Options() { addDoubleArrayOption("JST_SENSOR_COEFF", 2, false, jst_coeff); /*!\brief LAX_SENSOR_COEFF \n DESCRIPTION: 1st order artificial dissipation coefficients for the Lax-Friedrichs method. \ingroup Config*/ addDoubleOption("LAX_SENSOR_COEFF", Kappa_1st_Flow, 0.15); - /*!\brief USE_ACCURATE_FLUX_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for AUSM+up(2) and SLAU(2), exact Jacobians (for fixed wave speeds) for HLLC \ingroup Config*/ + /*!\brief USE_ACCURATE_FLUX_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for AUSM+up(2) and SLAU(2) \ingroup Config*/ addBoolOption("USE_ACCURATE_FLUX_JACOBIANS", Use_Accurate_Jacobians, false); + /*!\brief USE_EXACT_HLLC_JACOBIANS \n DESCRIPTION: Exact HLLC Jacobians with wave speeds held fixed \ingroup Config*/ + addBoolOption("USE_EXACT_HLLC_JACOBIANS", Use_Exact_HLLC_Jacobians, false); /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Improve the numerical properties (diagonal dominance) of the global Jacobian matrix, 3 to 4 is "optimum" (central schemes) \ingroup Config*/ addDoubleOption("CENTRAL_JACOBIAN_FIX_FACTOR", Cent_Jac_Fix_Factor, 4.0); /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Control numerical properties of the global Jacobian matrix using a multiplication factor for incompressible central schemes \ingroup Config*/ diff --git a/SU2_CFD/include/numerics/flow/convection/hllc.hpp b/SU2_CFD/include/numerics/flow/convection/hllc.hpp index 29f7ce727a34..8684862fa520 100644 --- a/SU2_CFD/include/numerics/flow/convection/hllc.hpp +++ b/SU2_CFD/include/numerics/flow/convection/hllc.hpp @@ -38,8 +38,9 @@ */ class CUpwHLLC_Flow final : public CNumerics { private: + friend struct HLLCJacobianTestAccess; bool implicit, dynamic_grid; - bool useAccurateJacobian; /*!< \brief Jacobians computed exactly for fixed wave speeds (USE_ACCURATE_FLUX_JACOBIANS). */ + bool useAccurateJacobian; /*!< \brief Jacobians computed exactly for fixed wave speeds (USE_EXACT_HLLC_JACOBIANS). */ unsigned short iDim, jDim, iVar, jVar; su2double *IntermediateState; @@ -81,19 +82,6 @@ class CUpwHLLC_Flow final : public CNumerics { * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. */ ResidualType<> ComputeResidual(const CConfig* config) override; - - /*! - * \brief Fix the left and right wave speeds (relative to the face), for the verification of the Jacobians, which - * treat them as constants. - * \param[in] active - Use the fixed values. - * \param[in] val_sL, val_sR - Left and right wave speeds. - */ - inline void SetFixedWaveSpeeds(bool active, su2double val_sL, su2double val_sR) { - fixedWaveSpeeds = active; - fixed_sL = val_sL; - fixed_sR = val_sR; - } - }; /*! @@ -105,8 +93,9 @@ class CUpwHLLC_Flow final : public CNumerics { */ class CUpwGeneralHLLC_Flow final : public CNumerics { private: + friend struct HLLCJacobianTestAccess; bool implicit, dynamic_grid; - bool useAccurateJacobian; /*!< \brief Jacobians computed exactly for fixed wave speeds (USE_ACCURATE_FLUX_JACOBIANS). */ + bool useAccurateJacobian; /*!< \brief Jacobians computed exactly for fixed wave speeds (USE_EXACT_HLLC_JACOBIANS). */ unsigned short iDim, jDim, iVar, jVar; su2double *IntermediateState; @@ -150,17 +139,6 @@ class CUpwGeneralHLLC_Flow final : public CNumerics { */ ResidualType<> ComputeResidual(const CConfig* config) override; - /*! - * \brief Fix the left and right wave speeds (relative to the face), for the verification of the Jacobians, which - * treat them as constants. - * \param[in] active - Use the fixed values. - * \param[in] val_sL, val_sR - Left and right wave speeds. - */ - inline void SetFixedWaveSpeeds(bool active, su2double val_sL, su2double val_sR) { - fixedWaveSpeeds = active; - fixed_sL = val_sL; - fixed_sR = val_sR; - } /*! * \brief Compute the Average quantities for a general fluid flux between two nodes i and j. diff --git a/SU2_CFD/src/numerics/flow/convection/hllc.cpp b/SU2_CFD/src/numerics/flow/convection/hllc.cpp index 2e315d6a1d91..131ee5b65fa7 100644 --- a/SU2_CFD/src/numerics/flow/convection/hllc.cpp +++ b/SU2_CFD/src/numerics/flow/convection/hllc.cpp @@ -80,7 +80,7 @@ CUpwHLLC_Flow::CUpwHLLC_Flow(unsigned short val_nDim, unsigned short val_nVar, c /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); - useAccurateJacobian = config->GetUse_Accurate_Jacobians(); + useAccurateJacobian = config->GetUse_Exact_HLLC_Jacobians(); Gamma = config->GetGamma(); @@ -292,13 +292,11 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Jacobians for the wave speeds held fixed (in the frame of the face). p* = rho_j (u_j - sR)(u_j - sM) + p_j = rho_i (u_i - sL)(u_i - sM) + p_i, so its derivative contains rho_j (sR - u_j) dsM/dU (or rho_i (sL - u_i) dsM/dU). By default the density of the other state is kept, as in the previous implementation (an inexact transcription - of the reference, see #2198). USE_ACCURATE_FLUX_JACOBIANS= YES gives the exact derivative. ---*/ + of the reference, see #2198). USE_EXACT_HLLC_JACOBIANS= YES gives the exact derivative. ---*/ const su2double rhoPStar_R = useAccurateJacobian ? Density_j : Density_i; const su2double rhoPStar_L = useAccurateJacobian ? Density_i : Density_j; - - if (sM > 0.0) { if (sL > 0.0) { @@ -597,7 +595,7 @@ CUpwGeneralHLLC_Flow::CUpwGeneralHLLC_Flow(unsigned short val_nDim, unsigned sho /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); - useAccurateJacobian = config->GetUse_Accurate_Jacobians(); + useAccurateJacobian = config->GetUse_Exact_HLLC_Jacobians(); Gamma = config->GetGamma(); @@ -837,13 +835,11 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Jacobians for the wave speeds held fixed (in the frame of the face). p* = rho_j (u_j - sR)(u_j - sM) + p_j = rho_i (u_i - sL)(u_i - sM) + p_i, so its derivative contains rho_j (sR - u_j) dsM/dU (or rho_i (sL - u_i) dsM/dU). By default the density of the other state is kept, as in the previous implementation (an inexact transcription - of the reference, see #2198). USE_ACCURATE_FLUX_JACOBIANS= YES gives the exact derivative. ---*/ + of the reference, see #2198). USE_EXACT_HLLC_JACOBIANS= YES gives the exact derivative. ---*/ const su2double rhoPStar_R = useAccurateJacobian ? Density_j : Density_i; const su2double rhoPStar_L = useAccurateJacobian ? Density_i : Density_j; - - if (sM > 0.0) { if (sL > 0.0) { diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp index 20a34724fadc..2dfe8b65cf70 100644 --- a/UnitTests/Common/CConfig_tests.cpp +++ b/UnitTests/Common/CConfig_tests.cpp @@ -101,3 +101,19 @@ TEST_CASE("Excessive multigrid levels are rejected", "[.MultigridLimitError]") { CConfig config(options, SU2_COMPONENT::SU2_CFD, false); FAIL("An unsupported multigrid hierarchy was accepted."); } + +TEST_CASE("HLLC and AUSM SLAU Jacobian options are independent", "[Config][HLLCOptions]") { + for (const auto accurate : {false, true}) { + for (const auto exactHLLC : {false, true}) { + std::stringstream options(base_options + "USE_ACCURATE_FLUX_JACOBIANS= " + (accurate ? "YES" : "NO") + + "\nUSE_EXACT_HLLC_JACOBIANS= " + (exactHLLC ? "YES" : "NO") + "\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + CHECK(config.GetUse_Accurate_Jacobians() == accurate); + CHECK(config.GetUse_Exact_HLLC_Jacobians() == exactHLLC); + } + } + std::stringstream options(base_options); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + CHECK_FALSE(config.GetUse_Accurate_Jacobians()); + CHECK_FALSE(config.GetUse_Exact_HLLC_Jacobians()); +} diff --git a/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp b/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp index b03edb3f1414..fc36c6cc9eba 100644 --- a/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp +++ b/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp @@ -1,7 +1,7 @@ /*! * \file hllc_jacobian_tests.cpp * \brief Finite-difference verification of the HLLC Jacobians (ideal and general gas, 2D and 3D), on fixed and moving - * faces, exact (USE_ACCURATE_FLUX_JACOBIANS) or default. + * faces, exact (USE_EXACT_HLLC_JACOBIANS) or default. * \version 8.5.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -32,6 +32,15 @@ #include #include "../../../SU2_CFD/include/numerics/flow/convection/hllc.hpp" +struct HLLCJacobianTestAccess { + template + static void SetFixedWaveSpeeds(Numerics& numerics, su2double sL, su2double sR) { + numerics.fixedWaveSpeeds = true; + numerics.fixed_sL = sL; + numerics.fixed_sR = sR; + } +}; + namespace { constexpr su2double gamma = 1.4; @@ -142,12 +151,12 @@ su2double ContactSpeed(unsigned short nDim, const State& Ui, const State& Uj, su std::unique_ptr MakeConfig(bool moving, bool accurate) { std::stringstream options; options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= EULER_IMPLICIT\n"; - options << "USE_ACCURATE_FLUX_JACOBIANS= " << (accurate ? "YES" : "NO") << "\n"; + options << "USE_EXACT_HLLC_JACOBIANS= " << (accurate ? "YES" : "NO") << "\n"; if (moving) options << "GRID_MOVEMENT= STEADY_TRANSLATION\nTRANSLATION_RATE= 0.0 0.0 0.0\n"; return std::make_unique(options, SU2_COMPONENT::SU2_CFD, false); } -/*--- With USE_ACCURATE_FLUX_JACOBIANS= YES the Jacobians are exact for fixed wave speeds in all branches. By default +/*--- With USE_EXACT_HLLC_JACOBIANS= YES the Jacobians are exact for fixed wave speeds in all branches. By default the derivative of p* is approximated (star states), the supersonic branches remain exact. ---*/ template void CheckAllBranches(unsigned short nDim, bool moving, su2double k, bool accurate) { @@ -186,7 +195,7 @@ void CheckAllBranches(unsigned short nDim, bool moving, su2double k, bool accura const su2double supersonicSpeeds[2] = {dir > 0 ? 0.02 : -2.0, dir > 0 ? 2.0 : -0.02}; const su2double starSpeeds[2] = {-2.0, 2.0}; for (const auto* s : {supersonicSpeeds, starSpeeds}) { - numerics.SetFixedWaveSpeeds(true, s[0], s[1]); + HLLCJacobianTestAccess::SetFixedWaveSpeeds(numerics, s[0], s[1]); const bool supersonic = (s[0] > 0.0) || (s[1] < 0.0); const su2double sM = ContactSpeed(nDim, Ui, Uj, k, 0.8 * k, n, w, s[0], s[1]); CAPTURE(nDim, moving, k, accurate, dir, s[0], s[1], sM); @@ -235,7 +244,7 @@ TEST_CASE("HLLC supersonic Jacobian uses reconstructed face enthalpy", "[HLLC-fa } numerics.SetNormal(normal); numerics.SetGridVel(gridVelocity, gridVelocity); - numerics.SetFixedWaveSpeeds(true, direction > 0 ? 2.0 : -4.0, direction > 0 ? 4.0 : -2.0); + HLLCJacobianTestAccess::SetFixedWaveSpeeds(numerics, direction > 0 ? 2.0 : -4.0, direction > 0 ? 4.0 : -2.0); const State state = Conservative(nDim, 1.2, velocity, 1.1, 0.0); for (su2double pressureOffset : {0.0, 0.1}) { auto evaluate = [&](const State& Ui, const State& Uj) { diff --git a/config_template.cfg b/config_template.cfg index 0d58bc62727a..8795186ab3d9 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1856,9 +1856,12 @@ LOW_MACH_PREC= NO % % Use numerically computed Jacobians for AUSM+up(2) and SLAU(2) % Slower per iteration but potentialy more stable and capable of higher CFL -% For HLLC: exact Jacobians for fixed wave speeds (the default approximation is more robust with inexact linear solves) USE_ACCURATE_FLUX_JACOBIANS= NO % +% Use exact HLLC Jacobians with wave speeds held fixed (NO, YES). +% The default retains the previous approximation; wave-speed derivatives are omitted in both modes. +USE_EXACT_HLLC_JACOBIANS= NO +% % Use the vectorized version of the selected numerical method (available for JST family and Roe). % SU2 should be compiled for an AVX or AVX512 architecture for best performance. % NOTE: Currently vectorization always used for schemes that support it. From 815aa5c610e5d0b6bb10086b4cf8a62a327ba5d1 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 22:28:45 +0200 Subject: [PATCH 12/12] Use deduced types for Roe entropy wave strengths --- SU2_CFD/src/numerics/flow/convection/roe.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/convection/roe.cpp b/SU2_CFD/src/numerics/flow/convection/roe.cpp index 37c697515cfb..3fd23d644584 100644 --- a/SU2_CFD/src/numerics/flow/convection/roe.cpp +++ b/SU2_CFD/src/numerics/flow/convection/roe.cpp @@ -336,7 +336,7 @@ void CUpwL2Roe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J } else { /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ - const su2double delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; + const auto delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; delta_wave[0] = UnitNormal[0]*delta_entropy + (delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1])*zeta; delta_wave[1] = UnitNormal[1]*delta_entropy + (delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2])*zeta; delta_wave[2] = UnitNormal[2]*delta_entropy + (delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0])*zeta; @@ -412,7 +412,7 @@ void CUpwLMRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J } else { /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ - const su2double delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; + const auto delta_entropy = delta_rho - delta_p/RoeSoundSpeed2; delta_wave[0] = UnitNormal[0]*delta_entropy + (delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1]); delta_wave[1] = UnitNormal[1]*delta_entropy + (delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2]); delta_wave[2] = UnitNormal[2]*delta_entropy + (delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0]); @@ -909,7 +909,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co } else { /*--- The first three columns of the 3D P matrix combine the entropy wave (weighted by the normal) with the shear waves, so their strengths are n_j (Delta rho - Delta p / c^2) + (Delta u x n)_j. ---*/ - const su2double delta_entropy = delta_rho - delta_p/(RoeSoundSpeed*RoeSoundSpeed); + const auto delta_entropy = delta_rho - delta_p/(RoeSoundSpeed*RoeSoundSpeed); delta_wave[0] = UnitNormal[0]*delta_entropy + delta_vel[1]*UnitNormal[2]-delta_vel[2]*UnitNormal[1]; delta_wave[1] = UnitNormal[1]*delta_entropy + delta_vel[2]*UnitNormal[0]-delta_vel[0]*UnitNormal[2]; delta_wave[2] = UnitNormal[2]*delta_entropy + delta_vel[0]*UnitNormal[1]-delta_vel[1]*UnitNormal[0];