From 77c905a32a15ab3b782f93f8668f4e2a6741d4c7 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Fri, 20 Jun 2025 11:28:25 +0100 Subject: [PATCH 1/8] sst under relaxation --- SU2_CFD/include/solvers/CTurbSSTSolver.hpp | 8 +++++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 35 ++++++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp index 5bdc3549aed3..352270f6011a 100644 --- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp @@ -53,6 +53,14 @@ class CTurbSSTSolver final : public CTurbSolver { CSolver **solver_container, const CConfig *config, unsigned short val_marker); + + /*! + * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over + * a nonlinear iteration for stability. + * \param[in] config - Definition of the particular problem. + */ + void ComputeUnderRelaxationFactor(const CConfig *config); + public: /*! * \brief Constructor. diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 9fee88d18760..5dc02de36232 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -1019,3 +1019,38 @@ void CTurbSSTSolver::SetUniformInlet(const CConfig* config, unsigned short iMark } } + +void CTurbSSTSolver::ComputeUnderRelaxationFactor(const CConfig *config) { + + const su2double allowableRatio = config->GetUnderRelax_SST(); + + SU2_OMP_FOR_STAT(omp_chunk_size) + /* Loop over the solution update given by relaxing the linear + system for this nonlinear iteration. */ + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + su2double localUnderRelaxation = 1.0; + + for (unsigned short iVar =0; iVar < nVar; iVar++) { + const unsigned long index = iPoint * nVar + iVar; + su2double ratio = fabs(LinSysSol[index])/(fabs(nodes->GetSolution(iPoint, iVar)) + EPS); + /* We impose a limit on the maximum percentage that the + turbulence variables can change over a nonlinear iteration. */ + if (ratio > allowableRatio) { + localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation); + } + + } + + /* Threshold the relaxation factor in the event that there is + a very small value. This helps avoid catastrophic crashes due + to non-realizable states by canceling the update. */ + + if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0; + + /* Store the under-relaxation factor for this point. */ + + nodes->SetUnderRelaxation(iPoint, localUnderRelaxation); + } + + END_SU2_OMP_FOR +} \ No newline at end of file From a210b9cc23d2b9170c8029b138efd8b871f26636 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Fri, 20 Jun 2025 11:29:04 +0100 Subject: [PATCH 2/8] add config options for UR values --- Common/include/CConfig.hpp | 21 +++++++++++++++++++++ Common/src/CConfig.cpp | 6 ++++++ 2 files changed, 27 insertions(+) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index af7aa99033b5..e5d96f3b689d 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -635,6 +635,9 @@ class CConfig { su2double Relaxation_Factor_Adjoint; /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */ su2double Relaxation_Factor_CHT; /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */ su2double EntropyFix_Coeff; /*!< \brief Entropy fix coefficient. */ + su2double SST_UnderRelaxation_Factor; /*!< \brief UnderRelaxation Factor for SST Turbulent Variables*/ + su2double SA_UnderRelaxation_Factor; /*!< \brief UnderRelaxation Factor for SST Turbulent Variables*/ + su2double Flow_UnderRelaxation_Factor; /*!< \brief UnderRelaxation Factor for SST Turbulent Variables*/ unsigned short nLocationStations, /*!< \brief Number of section cuts to make when outputting mesh and cp . */ nWingStations; /*!< \brief Number of section cuts to make when calculating internal volume. */ su2double Kappa_1st_AdjFlow, /*!< \brief Lax 1st order dissipation coefficient for adjoint flow equations (coarse multigrid levels). */ @@ -4207,6 +4210,24 @@ class CConfig { */ su2double GetRelaxation_Factor_CHT(void) const { return Relaxation_Factor_CHT; } + /*! + * \brief Get the under-relaxation for flow variables density and energy. + * \return under-relaxation for flow variables. + */ + su2double GetUnderRelax_Flow(void) const { return Flow_UnderRelaxation_Factor; } + + /*! + * \brief Get the under-relaxation for SA variable, nu_tilde. + * \return under-relaxation for SA variables. + */ + su2double GetUnderRelax_SA(void) const { return SA_UnderRelaxation_Factor; } + + /*! + * \brief Get the under-relaxation for SST turbulence variables k and omega. + * \return under-relaxation for SST variables. + */ + su2double GetUnderRelax_SST(void) const { return SST_UnderRelaxation_Factor; } + /*! * \brief Get the number of samples used in quasi-Newton methods. */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 5fb7220e94e4..0f6dadbfa0a0 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1887,6 +1887,12 @@ void CConfig::SetConfig_Options() { addEnumOption("DISCADJ_LIN_PREC", Kind_DiscAdj_Linear_Prec, Linear_Solver_Prec_Map, ILU); /* DESCRIPTION: Linear solver for the discete adjoint systems */ + addDoubleOption("UR_FACTOR_FLOW", Flow_UnderRelaxation_Factor, 0.2); + /* DESCRIPTION: Under-Relaxation Factor for Density and Energy Variables */ + addDoubleOption("UR_FACTOR_SA", SA_UnderRelaxation_Factor, 0.99); + /* DESCRIPTION: Under-Relaxation Factor for SA Turbulence Variables */ + addDoubleOption("UR_FACTOR_SST", SST_UnderRelaxation_Factor, 1.0); + /* DESCRIPTION: Under-Relaxation Factor for SST Turbulence Variables */ /*!\par CONFIG_CATEGORY: Convergence\ingroup Config*/ /*--- Options related to convergence ---*/ From 5c33dd08c1149e1fa4eed9be97e8738c44bd6367 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Fri, 20 Jun 2025 11:30:14 +0100 Subject: [PATCH 3/8] ur value based on config option (default is same as before) --- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index d3673e1b763a..39b35f7d1c4e 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -558,7 +558,7 @@ void CFVMFlowSolverBase::ComputeUnderRelaxationFactor(const CConfig* confi /* Loop over the solution update given by relaxing the linear system for this nonlinear iteration. */ - const su2double allowableRatio = 0.2; + const su2double allowableRatio = config->GetUnderRelax_Flow(); SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { From 3420163c51dcf5711f1c8888b61d0e3fed0a1557 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Fri, 20 Jun 2025 11:30:36 +0100 Subject: [PATCH 4/8] ur value based on config option (default is same as before) --- SU2_CFD/src/solvers/CTurbSASolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 4a0d54a95176..a0f3bfe234d4 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -1560,7 +1560,7 @@ void CTurbSASolver::ComputeUnderRelaxationFactor(const CConfig *config) { system for this nonlinear iteration. */ su2double localUnderRelaxation = 1.00; - const su2double allowableRatio = 0.99; + const su2double allowableRatio = config->GetUnderRelax_SA(); SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { From fb6bd4cbe74e70f7c6fdd443c0dd8ad4d71cc7f2 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Wed, 25 Jun 2025 12:55:52 +0100 Subject: [PATCH 5/8] stop CD residual accumulation as source --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 4fa9080e1ad5..37e53807628a 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -924,7 +924,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Cross diffusion ---*/ - Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; + // Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; /*--- Contribution due to 2D axisymmetric formulation ---*/ From e60372c2205f7e6cada812945414d3d9b2ac245e Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Wed, 25 Jun 2025 12:56:19 +0100 Subject: [PATCH 6/8] option 1 with frozen coefficient jacobians --- .../numerics/turbulent/turb_diffusion.hpp | 31 ++++++++++++++----- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp index b8edcd9007a7..6af1c6570949 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp @@ -223,28 +223,43 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { const su2double sigma_kine_j = F1_j*sigma_k1 + (1.0 - F1_j)*sigma_k2; const su2double sigma_omega_i = F1_i*sigma_om1 + (1.0 - F1_i)*sigma_om2; const su2double sigma_omega_j = F1_j*sigma_om1 + (1.0 - F1_j)*sigma_om2; + const su2double lambda_i = 2 * (1 - F1_i) * Density_i * sigma_omega_i / ScalarVar_i[1]; + const su2double lambda_j = 2 * (1 - F1_j) * Density_j * sigma_omega_j / ScalarVar_j[1]; + const su2double lambda_ij = 0.5 * (lambda_i + lambda_j); /*--- Compute mean effective dynamic viscosity ---*/ const su2double diff_i_kine = Laminar_Viscosity_i + sigma_kine_i*Eddy_Viscosity_i; const su2double diff_j_kine = Laminar_Viscosity_j + sigma_kine_j*Eddy_Viscosity_j; - const su2double diff_i_omega = Laminar_Viscosity_i + sigma_omega_i*Eddy_Viscosity_i; - const su2double diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j; + const su2double diff_i_omega = Laminar_Viscosity_i + sigma_omega_i*Eddy_Viscosity_i + lambda_i * ScalarVar_i[0]; + const su2double diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j + lambda_j * ScalarVar_j[0]; + + //non-conservative term: + const su2double diff_omega_t1 = 0.5 * (diff_i_omega + diff_j_omega); + const su2double diff_omega_t2 = -ScalarVar_i[0] * lambda_ij; const su2double diff_kine = 0.5*(diff_i_kine + diff_j_kine); - const su2double diff_omega = 0.5*(diff_i_omega + diff_j_omega); + const su2double diff_omega = diff_omega_t1 + diff_omega_t2; Flux[0] = diff_kine*Proj_Mean_GradScalarVar[0]; Flux[1] = diff_omega*Proj_Mean_GradScalarVar[1]; - /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ + /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients + * Using chain rule for dFlux[1]/dx_{i.j} = d(diff_coef)/dx_{i,j}*Proj + d(Proj)/dx_{i,j}*diff_ceof---*/ + if (implicit) { const su2double proj_on_rho_i = proj_vector_ij/Density_i; - Jacobian_i[0][0] = -diff_kine*proj_on_rho_i; Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_omega*proj_on_rho_i; + const su2double dlambda_domega_i = - (2 * (1 - F1_i) * Density_i * (sigma_omega_i))/pow(ScalarVar_i[1], 2); + Jacobian_i[0][0] = -diff_kine*proj_on_rho_i; + Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; //dF_w/dk_i + Jacobian_i[1][1] = (diff_omega) * -proj_on_rho_i; //lambda not frozen const su2double proj_on_rho_j = proj_vector_ij/Density_j; - Jacobian_j[0][0] = diff_kine*proj_on_rho_j; Jacobian_j[0][1] = 0.0; - Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_omega*proj_on_rho_j; + const su2double dlambda_domega_j = - (2 * (1 - F1_j) * Density_j * (sigma_omega_j))/pow(ScalarVar_j[1], 2); + Jacobian_j[0][0] = diff_kine*proj_on_rho_j; + Jacobian_j[0][1] = 0.0; + Jacobian_j[1][0] = 0.0; + Jacobian_j[1][1] = (diff_omega) * proj_on_rho_j; //frozen diff_ceoff } } From 6ecadd59c9757d06c14752b17b0335a16b207fc0 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Wed, 25 Jun 2025 12:56:51 +0100 Subject: [PATCH 7/8] update scalar solver --- SU2_CFD/include/solvers/CScalarSolver.hpp | 1 + SU2_CFD/include/solvers/CScalarSolver.inl | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index 0fc61aa7ab25..bbefefc5bf4e 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -77,6 +77,7 @@ class CScalarSolver : public CSolver { /*--- Edge fluxes for reducer strategy (see the notes in CEulerSolver.hpp). ---*/ CSysVector EdgeFluxes; /*!< \brief Flux across each edge. */ + CSysVector EdgeFluxes_Diff; /*!< \brief Flux difference across ij and ji for non-conservative discretisation. */ /*! * \brief The highest level in the variable hierarchy this solver can safely use. diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 314354116ff1..f7029e1c9ad9 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -348,8 +348,12 @@ void CScalarSolver::SumEdgeFluxes(CGeometry* geometry) { for (auto iEdge : geometry->nodes->GetEdges(iPoint)) { if (iPoint == geometry->edges->GetNode(iEdge, 0)) LinSysRes.AddBlock(iPoint, EdgeFluxes.GetBlock(iEdge)); - else + else { LinSysRes.SubtractBlock(iPoint, EdgeFluxes.GetBlock(iEdge)); + if (EdgeFluxes_Diff.GetLocSize() > 0) { + LinSysRes.SubtractBlock(iPoint, EdgeFluxes_Diff.GetBlock(iEdge)); + } + } } } END_SU2_OMP_FOR From 786b404f4532e3255a8d8bfe84c01cf12a402696 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Wed, 25 Jun 2025 12:57:07 +0100 Subject: [PATCH 8/8] update SST solver --- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 93 +++++++++++++++++++++++++- 1 file changed, 91 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 5dc02de36232..b29b6622ddec 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -75,8 +75,10 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); System.SetxIsZero(true); - if (ReducerStrategy) + if (ReducerStrategy) { EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + EdgeFluxes_Diff.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + } /*--- Initialize the BGS residuals in multizone problems. ---*/ if (multizone){ @@ -295,8 +297,95 @@ void CTurbSSTSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry }; /*--- Now instantiate the generic implementation with the functor above. ---*/ + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + CFlowVariable* flowNodes = solver_container[FLOW_SOL] ? + su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()) : nullptr; + + /*--- Points in edge ---*/ + auto iPoint = geometry->edges->GetNode(iEdge, 0); + auto jPoint = geometry->edges->GetNode(iEdge, 1); + + /*--- Helper function to compute the flux ---*/ + auto ComputeFlux = [&](unsigned long point_i, unsigned long point_j, const su2double* normal) { + numerics->SetCoord(geometry->nodes->GetCoord(point_i),geometry->nodes->GetCoord(point_j)); + numerics->SetNormal(normal); + + if (flowNodes) { + numerics->SetPrimitive(flowNodes->GetPrimitive(point_i), flowNodes->GetPrimitive(point_j)); + } + + numerics->SetScalarVar(nodes->GetSolution(point_i), nodes->GetSolution(point_j)); + numerics->SetScalarVarGradient(nodes->GetGradient(point_i), nodes->GetGradient(point_j)); + + return numerics->ComputeResidual(config); + }; + + SolverSpecificNumerics(iPoint, jPoint); + + /*--- Compute fluxes and jacobians i->j ---*/ + const su2double* normal = geometry->edges->GetNormal(iEdge); + auto residual_ij = ComputeFlux(iPoint, jPoint, normal); + if (ReducerStrategy) { + EdgeFluxes.SubtractBlock(iEdge, residual_ij); + EdgeFluxes_Diff.SetBlock(iEdge, residual_ij); + if (implicit) { + auto* Block_ij = Jacobian.GetBlock(iPoint, jPoint); + auto* Block_ji = Jacobian.GetBlock(jPoint, iPoint); + + for (int iVar=0; iVari ---*/ + su2double flipped_normal[MAXNDIM]; + for (int iDim=0; iDim