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 ---*/ 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 } } 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 ---*/ 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++) { 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 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/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++) { diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 9fee88d18760..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; iDimGetUnderRelax_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