From 876c5b26e375409a804b8e28b0d62fe613305fe0 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 5 Oct 2025 23:26:04 +0200 Subject: [PATCH] initial commit from Bellosta, Abergo and Nishikawa code --- Common/include/CConfig.hpp | 7 +++ Common/include/option_structure.hpp | 16 ++++++ Common/src/CConfig.cpp | 1 + SU2_CFD/include/numerics/CNumerics.hpp | 27 +++++++-- .../include/numerics/flow/flow_diffusion.hpp | 19 +++++-- SU2_CFD/include/numerics/heat.hpp | 2 +- .../numerics/scalar/scalar_diffusion.hpp | 7 ++- .../numerics/species/species_diffusion.hpp | 2 +- .../turbulent/transition/trans_diffusion.hpp | 2 +- .../numerics/turbulent/turb_diffusion.hpp | 12 ++-- .../numerics_simd/flow/diffusion/common.hpp | 19 ++++++- .../flow/diffusion/viscous_fluxes.hpp | 30 ++++++++-- .../include/numerics_simd/flow/variables.hpp | 30 ++++++++++ SU2_CFD/src/drivers/CDriver.cpp | 42 +++++++------- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 56 ++++++++++++++++--- SU2_CFD/src/numerics/radiation.cpp | 2 +- 16 files changed, 215 insertions(+), 59 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index abbb0f401da6..9cddc4d977a5 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -569,6 +569,7 @@ class CConfig { Kind_Centered_Template; /*!< \brief Centered scheme for the template model. */ + VISCOUS_GRAD_CORR Kind_Visc_Corr; /*!< \brief Gradient correction damping. */ FEM_SHOCK_CAPTURING_DG Kind_FEM_Shock_Capturing_DG; /*!< \brief Shock capturing method for the FEM DG solver. */ BGS_RELAXATION Kind_BGS_RelaxMethod; /*!< \brief Kind of relaxation method for Block Gauss Seidel method in FSI problems. */ bool ReconstructionGradientRequired; /*!< \brief Enable or disable a second gradient calculation for upwind reconstruction only. */ @@ -4595,6 +4596,12 @@ class CConfig { */ UPWIND GetKind_Upwind(void) const { return Kind_Upwind; } + /*! + * \brief Get kind of gradient correction damping term for the viscous terms. + * \return Kind of gradient correction damping term. + */ + VISCOUS_GRAD_CORR GetKind_ViscousGradCorr(void) const { return Kind_Visc_Corr; } + /*! * \brief Get if the upwind scheme used MUSCL or not. * \note This is the information that the code will use, the method will diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bf99257b1387..c50df312a093 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -878,6 +878,22 @@ static const MapType Upwind_Map = { MakePair("LAX-FRIEDRICH", UPWIND::LAX_FRIEDRICH) }; +/*! + * \brief Types of viscous gradient correction spatial discretizations + */ +enum class VISCOUS_GRAD_CORR { + NONE, /*!< \brief No correction is used. */ + EDGE_NORMAL, /*!< \brief Edge Normal Correction. */ + FACE_TANGENT, /*!< \brief Face Tangent correction. */ + ALPHA_DAMPING /*!< \brief Alpha Damping with a=5/4. */ +}; +static const MapType Viscous_Grad_Corr_Map = { + MakePair("NONE", VISCOUS_GRAD_CORR::NONE) + MakePair("EDGE_NORMAL", VISCOUS_GRAD_CORR::EDGE_NORMAL) + MakePair("FACE_TANGENT", VISCOUS_GRAD_CORR::FACE_TANGENT) + MakePair("ALPHA_DAMPING", VISCOUS_GRAD_CORR::ALPHA_DAMPING) +}; + /*! * \brief Types of FEM spatial discretizations */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index a2a529917bf3..f647a1339f7a 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1693,6 +1693,7 @@ void CConfig::SetConfig_Options() { addStringDoubleListOption("MARKER_HEATFLUX", nMarker_HeatFlux, Marker_HeatFlux, Heat_Flux); /*!\brief INTEGRATED_HEATFLUX \n DESCRIPTION: Prescribe Heatflux in [W] instead of [W/m^2] \ingroup Config \default false */ addBoolOption("INTEGRATED_HEATFLUX", Integrated_HeatFlux, false); + addEnumOption("VISC_GRAD_CORR", Kind_Visc_Corr, Viscous_Grad_Corr_Map, VISCOUS_GRAD_CORR::EDGE_NORMAL); /*!\brief MARKER_HEATTRANSFER DESCRIPTION: Heat flux with specified heat transfer coefficient boundary marker(s)\n * Format: ( Heat transfer marker, heat transfer coefficient, wall temperature (static), ... ) \ingroup Config */ addExhaustOption("MARKER_HEATTRANSFER", nMarker_HeatTransfer, Marker_HeatTransfer, HeatTransfer_Coeff, HeatTransfer_WallTemp); diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index d5fb6a0ec11e..7b4c932d6013 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -654,20 +654,39 @@ class CNumerics { FORCEINLINE static su2double ComputeProjectedGradient(int nDim, int nVar, const Vec1& normal, const Vec1& coord_i, const Vec1& coord_j, const Mat& grad_i, const Mat& grad_j, - bool correct, + const VISCOUS_GRAD_CORR correct, const Vec2& var_i, const Vec2& var_j, su2double* projNormal, su2double* projCorrected) { assert(nDim == 2 || nDim == 3); nDim = (nDim > 2)? 3 : 2; su2double edgeVec[MAXNDIM], dist_ij_2 = 0.0, proj_vector_ij = 0.0; + su2double diss = 0.0; + su2double Area2 = 0.0; + for (int iDim = 0; iDim < nDim; ++iDim) { + Area2 += pow(normal[iDim],2); + } for (int iDim = 0; iDim < nDim; iDim++) { edgeVec[iDim] = coord_j[iDim] - coord_i[iDim]; dist_ij_2 += pow(edgeVec[iDim], 2); proj_vector_ij += edgeVec[iDim] * normal[iDim]; } - proj_vector_ij /= max(dist_ij_2,EPS); + + switch (correct) { + case VISCOUS_GRAD_CORR::EDGE_NORMAL: + diss = proj_vector_ij / max(dist_ij_2,EPS*EPS); + break; + case VISCOUS_GRAD_CORR::FACE_TANGENT: + diss = Area2 / max(proj_vector_ij,EPS); + break; + case VISCOUS_GRAD_CORR::ALPHA_DAMPING: + const su2double alpha = 4.0 / 3.0; + diss = alpha * Area2 / max(abs(proj_vector_ij),EPS); + break; + } + + //proj_vector_ij /= max(dist_ij_2,EPS); /*--- Mean gradient approximation. ---*/ for (int iVar = 0; iVar < nVar; iVar++) { @@ -677,11 +696,11 @@ class CNumerics { for (int iDim = 0; iDim < nDim; iDim++) { su2double meanGrad = 0.5 * (grad_i[iVar][iDim] + grad_j[iVar][iDim]); projNormal[iVar] += meanGrad * normal[iDim]; - if (correct) edgeProj += meanGrad * edgeVec[iDim]; + edgeProj += meanGrad * edgeVec[iDim]; } projCorrected[iVar] = projNormal[iVar]; - if (correct) projCorrected[iVar] -= (edgeProj - (var_j[iVar]-var_i[iVar])) * proj_vector_ij; + projCorrected[iVar] += ((var_j[iVar]-var_i[iVar]) - edgeProj) * diss; } return proj_vector_ij; diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp index 59cced23e6bc..b1cf7e4fd9ec 100644 --- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp +++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp @@ -44,7 +44,7 @@ class CAvgGrad_Base : public CNumerics { protected: const unsigned short nPrimVar; /*!< \brief The size of the primitive variable array used in the numerics class. */ - const bool correct_gradient; /*!< \brief Apply a correction to the gradient term */ + const VISCOUS_GRAD_CORR correct_gradient; /*!< \brief Apply a correction to the gradient term */ bool implicit = false; /*!< \brief Implicit calculus. */ su2double heat_flux_vector[MAXNDIM] = {0.0}, /*!< \brief Flux of total energy due to molecular and turbulent diffusion */ @@ -164,6 +164,15 @@ class CAvgGrad_Base : public CNumerics { su2double val_dist_ij_2, const unsigned short val_nPrimVar); + void GradientCorrection(su2double** GradPrimVar, + const su2double* val_PrimVar_i, + const su2double* val_PrimVar_j, + const su2double* val_edge_vector, + su2double val_dist_ij_2, + const su2double* diss, + const unsigned short val_nPrimVar); + + public: /*! @@ -176,7 +185,7 @@ class CAvgGrad_Base : public CNumerics { */ CAvgGrad_Base(unsigned short val_nDim, unsigned short val_nVar, unsigned short val_nPrimVar, - bool val_correct_grad, const CConfig* config); + VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config); /*! * \brief Destructor of the class. @@ -251,7 +260,7 @@ class CAvgGrad_Flow final : public CAvgGrad_Base { * \param[in] config - Definition of the particular problem. */ CAvgGrad_Flow(unsigned short val_nDim, unsigned short val_nVar, - bool val_correct_grad, const CConfig* config); + VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config); /*! * \brief Compute the viscous flow residual using an average of gradients. @@ -330,7 +339,7 @@ class CAvgGradInc_Flow final : public CAvgGrad_Base { * \param[in] config - Definition of the particular problem. */ CAvgGradInc_Flow(unsigned short val_nDim, unsigned short val_nVar, - bool val_correct_grad, const CConfig* config); + VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config); /*! * \brief Compute the viscous flow residual using an average of gradients. @@ -383,7 +392,7 @@ class CGeneralAvgGrad_Flow final : public CAvgGrad_Base { * \param[in] val_correct_grad - Apply a correction to the gradient * \param[in] config - Definition of the particular problem. */ - CGeneralAvgGrad_Flow(unsigned short val_nDim, unsigned short val_nVar, bool val_correct_grad, const CConfig* config); + CGeneralAvgGrad_Flow(unsigned short val_nDim, unsigned short val_nVar, VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config); /*! * \brief Compute the viscous flow residual using an average of gradients. diff --git a/SU2_CFD/include/numerics/heat.hpp b/SU2_CFD/include/numerics/heat.hpp index dceec7c510b9..63ead3f43659 100644 --- a/SU2_CFD/include/numerics/heat.hpp +++ b/SU2_CFD/include/numerics/heat.hpp @@ -80,7 +80,7 @@ class CAvgGrad_Heat final : public CAvgGrad_Scalar { * \param[in] config - Definition of the particular problem. * \param[in] correct - Whether to correct the gradient. */ - CAvgGrad_Heat(unsigned short val_nDim, const CConfig *config, bool correct) + CAvgGrad_Heat(unsigned short val_nDim, const CConfig *config, VISCOUS_GRAD_CORR correct) : CAvgGrad_Scalar(val_nDim, 1, correct, config) {} private: diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp index 13ee01f7eda6..48115e0f664a 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp @@ -69,7 +69,8 @@ class CAvgGrad_Scalar : public CNumerics { su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */ su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */ - const bool correct_gradient = false, incompressible = false; + const bool incompressible = false; + const VISCOUS_GRAD_CORR correct_gradient; /*! * \brief A pure virtual function; Adds any extra variables to AD @@ -91,7 +92,7 @@ class CAvgGrad_Scalar : public CNumerics { * \param[in] correct_gradient - Whether to correct gradient for skewness. * \param[in] config - Definition of the particular problem. */ - CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, + CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, VISCOUS_GRAD_CORR correct_grad, const CConfig* config) : CNumerics(val_nDim, val_nVar, config), idx(val_nDim, config->GetnSpecies()), @@ -123,7 +124,7 @@ class CAvgGrad_Scalar : public CNumerics { AD::SetPreaccIn(Normal, nDim); AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); AD::SetPreaccIn(ScalarVar_Grad_j, nVar, nDim); - if (correct_gradient) { + if (correct_gradient != VISCOUS_GRAD_CORR::NONE) { AD::SetPreaccIn(ScalarVar_i, nVar); AD::SetPreaccIn(ScalarVar_j, nVar); } diff --git a/SU2_CFD/include/numerics/species/species_diffusion.hpp b/SU2_CFD/include/numerics/species/species_diffusion.hpp index 3200f5125313..689cdfe8784e 100644 --- a/SU2_CFD/include/numerics/species/species_diffusion.hpp +++ b/SU2_CFD/include/numerics/species/species_diffusion.hpp @@ -104,7 +104,7 @@ class CAvgGrad_Species final : public CAvgGrad_Scalar { * \param[in] correct_grad - Whether to correct gradient for skewness. * \param[in] config - Definition of the particular problem. */ - CAvgGrad_Species(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config) + CAvgGrad_Species(unsigned short val_nDim, unsigned short val_nVar, VISCOUS_GRAD_CORR correct_grad, const CConfig* config) : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), turbulence(config->GetKind_Turb_Model() != TURB_MODEL::NONE) {} }; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp index 3c4ce9b6faee..795cb80dac18 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp @@ -98,7 +98,7 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar { * \param[in] correct_grad - Whether to correct gradient for skewness. * \param[in] config - Definition of the particular problem. */ - CAvgGrad_TransLM(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config) + CAvgGrad_TransLM(unsigned short val_nDim, unsigned short val_nVar, VISCOUS_GRAD_CORR correct_grad, const CConfig* config) : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config){ } diff --git a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp index 5d58f76b671c..3e8ad1840e62 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp @@ -115,7 +115,7 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { * \param[in] config - Definition of the particular problem. */ CAvgGrad_TurbSA(unsigned short val_nDim, unsigned short val_nVar, - bool correct_grad, const CConfig* config) + VISCOUS_GRAD_CORR correct_grad, const CConfig* config) : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) {} }; @@ -208,7 +208,7 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { * \param[in] config - Definition of the particular problem. */ CAvgGrad_TurbSA_Neg(unsigned short val_nDim, unsigned short val_nVar, - bool correct_grad, const CConfig* config) + VISCOUS_GRAD_CORR correct_grad, const CConfig* config) : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config) {} }; @@ -275,7 +275,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { /*--- We aim to treat the cross-diffusion as a diffusion term rather than a source term. * Re-writing the cross-diffusion contribution as λ/w ∇w ∇k, where λ = (2 (1- F1) ρ σ_ω2) - * and expanding using the product rule for divergence theorem gives: ∇(w λ/w ∇k) - w ∇(λ/w ∇k). + * and expanding using the product rule for divergence theorem gives: ∇(w λ/w ∇k) - w ∇(λ/w ∇k). * Discretising using FVM, gives: (λ)_ij ∇k - w_c (λ/w)_ij ∇k. where w_c is the cell centre value ---*/ const su2double lambda_i = 2 * (1 - F1_i) * Density_i * sigma_omega_i; @@ -284,7 +284,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { const su2double w_ij = 0.5 * (ScalarVar_i[1] + ScalarVar_j[1]); const su2double diff_omega_T2 = lambda_ij; - + const su2double diff_omega_T3 = -ScalarVar_i[1] * lambda_ij/w_ij; Flux[0] = diff_kine*Proj_Mean_GradScalarVar[0]; @@ -314,7 +314,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { Jacobian_j[0][1] = 0.0; Jacobian_j[1][0] = (diff_omega_T2 + diff_omega_T3)*proj_on_rho_j; Jacobian_j[1][1] = proj_on_rho_j * diff_omega_T1 + 2*lambda_ij*ScalarVar_i[1]/pow(ScalarVar_i[1]+ScalarVar_j[1],2) * Proj_Mean_GradScalarVar[0]; - } + } } } @@ -328,7 +328,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { * \param[in] config - Definition of the particular problem. */ CAvgGrad_TurbSST(unsigned short val_nDim, unsigned short val_nVar, - const su2double* constants, bool correct_grad, const CConfig* config) + const su2double* constants, VISCOUS_GRAD_CORR correct_grad, const CConfig* config) : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), sigma_k1(constants[0]), sigma_k2(constants[1]), diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp index 885c7494e7f6..26106a0753b1 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp @@ -49,18 +49,31 @@ FORCEINLINE MatrixDbl averageGradient(Int iPoint, Int jPoint, return avgGrad; } +template +FORCEINLINE CCompressibleSecondary averageSecondary(Int iPoint, Int jPoint, + const SecondaryType& secondary) { + + CCompressibleSecondary out; + auto second_i = gatherVariables(iPoint, secondary); + auto second_j = gatherVariables(jPoint, secondary); + for (size_t iVar = 0; iVar < nSecVar; ++iVar) { + out.all(iVar) = 0.5 * (second_i(iVar) + second_j(iVar)); + } + return out; +} + /*! * \brief Correct average gradient with the directional derivative to avoid decoupling. */ template FORCEINLINE void correctGradient(const PrimitiveType& V, const VectorDbl& vector_ij, - Double dist2_ij, + const VectorDbl& diss, MatrixDbl& avgGrad) { for (size_t iVar = 0; iVar < nVar; ++iVar) { - Double corr = (dot(avgGrad[iVar],vector_ij) - V.j.all(iVar) + V.i.all(iVar)) / dist2_ij; + Double corr = ( V.j.all(iVar) - V.i.all(iVar) - dot(avgGrad[iVar],vector_ij)); for (size_t iDim = 0; iDim < nDim; ++iDim) { - avgGrad(iVar,iDim) -= corr * vector_ij(iDim); + avgGrad(iVar,iDim) += corr * diss(iDim); } } } diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp index be271f7fd1c2..38d5225a79f0 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -71,12 +71,16 @@ template class CCompressibleViscousFluxBase : public CNumericsSIMD { protected: static constexpr size_t nDim = NDIM; - static constexpr size_t nPrimVarGrad = nDim+1; + static constexpr size_t nPrimVarGrad = nDim+3; + static constexpr size_t nSeconVar = 8; const su2double gamma; const su2double gasConst; const su2double prandtlTurb; - const bool correct; + //const su2double cp; + const bool correct_EN; + const bool correct_FT; + const bool correct_AD; const bool useSA_QCR; const bool wallFun; const bool uq; @@ -96,7 +100,9 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { gamma(config.GetGamma()), gasConst(config.GetGas_ConstantND()), prandtlTurb(config.GetPrandtl_Turb()), - correct(iMesh == MESH_0), + correct_EN(iMesh == MESH_0 && config.GetKind_ViscousGradCorr() == VISCOUS_GRAD_CORR::EDGE_NORMAL), + correct_FT(iMesh == MESH_0 && config.GetKind_ViscousGradCorr() == VISCOUS_GRAD_CORR::FACE_TANGENT), + correct_AD(iMesh == MESH_0 && config.GetKind_ViscousGradCorr() == VISCOUS_GRAD_CORR::ALPHA_DAMPING), useSA_QCR(config.GetSAParsedOptions().qcr2000), wallFun(config.GetWall_Functions()), uq(config.GetSSTParsedOptions().uq), @@ -135,6 +141,7 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { const auto& solution = static_cast(solution_); const auto& gradient = solution.GetGradient_Primitive(); + const auto& secondary = solution.GetSecondary(); /*--- Compute distance and handle zero without "ifs" by making it large. ---*/ @@ -145,7 +152,22 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { /*--- Compute the corrected mean gradient. ---*/ auto avgGrad = averageGradient(iPoint, jPoint, gradient); - if(correct) correctGradient(V, vector_ij, dist2_ij, avgGrad); + auto avgSecond = averageSecondary(iPoint, jPoint, secondary); + + Double eDn = dot(vector_ij,unitNormal); + mask = eDn < EPS; + eDn += mask / (EPS); + const Double eDotN = Double(1.0) / eDn; + const Double alpha = Double(4.0) / 3.0; + + VectorDbl diss; + + if(correct_EN) for (int iDim = 0; iDim < nDim; ++iDim) diss(iDim) = vector_ij(iDim) / dist2_ij; + else if(correct_FT) for (int iDim = 0; iDim < nDim; ++iDim) diss(iDim) = unitNormal(iDim) * eDotN; + else if(correct_AD) for (int iDim = 0; iDim < nDim; ++iDim) diss(iDim) = unitNormal(iDim) * alpha * abs(eDotN); + + if (correct_EN || correct_FT || correct_AD) + correctGradient(V, vector_ij, diss, avgGrad); /*--- Stress and heat flux tensors. ---*/ diff --git a/SU2_CFD/include/numerics_simd/flow/variables.hpp b/SU2_CFD/include/numerics_simd/flow/variables.hpp index 12dfb560c23e..a9e04c84a56b 100644 --- a/SU2_CFD/include/numerics_simd/flow/variables.hpp +++ b/SU2_CFD/include/numerics_simd/flow/variables.hpp @@ -83,6 +83,36 @@ struct CCompressibleConservatives { FORCEINLINE const Double* momentum() const { return &momentum(0); } }; +/*! + * \brief Type to store compressible conservative (i.e. solution) variables. + */ +template +struct CCompressibleSecondary { + /*!< \brief Secondary variables (dPdrho_e, dPde_rho, dTdrho_e, dTde_rho, dmudrho_T, dmudT_rho, dktdrho_T, dktdT_rho) + * in compressible (Euler: 2, NS: 8) flows. */ + static constexpr size_t nVar = nVar_; + VectorDbl all; + + FORCEINLINE Double& dPdrho_e() { return all(0); } + FORCEINLINE Double& dPde_rho() { return all(1); } + FORCEINLINE Double& dTdrho_e() { return all(2); } + FORCEINLINE Double& dTde_rho() { return all(3); } + FORCEINLINE Double& dmudrho_T() { return all(4); } + FORCEINLINE Double& dmudT_rho() { return all(5); } + FORCEINLINE Double& dktdrho_T() { return all(6); } + FORCEINLINE Double& dktdT_rho() { return all(7); } + + FORCEINLINE const Double& dPdrho_e() const { return all(0); } + FORCEINLINE const Double& dPde_rho() const { return all(1); } + FORCEINLINE const Double& dTdrho_e() const { return all(2); } + FORCEINLINE const Double& dTde_rho() const { return all(3); } + FORCEINLINE const Double& dmudrho_T() const { return all(4); } + FORCEINLINE const Double& dmudT_rho() const { return all(5); } + FORCEINLINE const Double& dktdrho_T() const { return all(6); } + FORCEINLINE const Double& dktdT_rho() const { return all(7); } + +}; + /*! * \brief Primitive to conservative conversion. */ diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 7e17ed9c6015..a1bee09a8012 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1256,13 +1256,13 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { if (spalart_allmaras) { if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG) { - numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, true, config); + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, config->GetKind_ViscousGradCorr(), config); } else { - numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, true, config); + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, config->GetKind_ViscousGradCorr(), config); } } else if (menter_sst) - numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, true, config); + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, config->GetKind_ViscousGradCorr(), config); } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -1286,14 +1286,14 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG) { - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, VISCOUS_GRAD_CORR::NONE, config); } else { - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, VISCOUS_GRAD_CORR::NONE, config); } } else if (menter_sst) { numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, false, + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, VISCOUS_GRAD_CORR::NONE, config); } } @@ -1341,7 +1341,7 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse /*--- Definition of the viscous scheme for each equation and mesh level ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, true, config); + if (LM) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, config->GetKind_ViscousGradCorr(), config); } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -1359,7 +1359,7 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { if (LM) { numerics[iMGlevel][TRANS_SOL][conv_bound_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); - numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, false, config); + numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, VISCOUS_GRAD_CORR::NONE, config); } } } @@ -1404,8 +1404,8 @@ void CDriver::InstantiateSpeciesNumerics(unsigned short nVar_Species, int offset /*--- Definition of the viscous scheme for each equation and mesh level ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - numerics[iMGlevel][SPECIES_SOL][visc_term] = new CAvgGrad_Species(nDim, nVar_Species, true, config); - numerics[iMGlevel][SPECIES_SOL][visc_bound_term] = new CAvgGrad_Species(nDim, nVar_Species, false, config); + numerics[iMGlevel][SPECIES_SOL][visc_term] = new CAvgGrad_Species(nDim, nVar_Species, config->GetKind_ViscousGradCorr(), config); + numerics[iMGlevel][SPECIES_SOL][visc_bound_term] = new CAvgGrad_Species(nDim, nVar_Species, VISCOUS_GRAD_CORR::NONE, config); } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -1801,36 +1801,36 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver if (ideal_gas) { /*--- Compressible flow Ideal gas ---*/ - numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGrad_Flow(nDim, nVar_Flow, true, config); + numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGrad_Flow(nDim, nVar_Flow, config->GetKind_ViscousGradCorr(), config); for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGrad_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGrad_Flow(nDim, nVar_Flow, VISCOUS_GRAD_CORR::NONE, config); /*--- Definition of the boundary condition method ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGrad_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGrad_Flow(nDim, nVar_Flow, VISCOUS_GRAD_CORR::NONE, config); } else { /*--- Compressible flow Real gas ---*/ - numerics[MESH_0][FLOW_SOL][visc_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, true, config); + numerics[MESH_0][FLOW_SOL][visc_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, config->GetKind_ViscousGradCorr(), config); for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, VISCOUS_GRAD_CORR::NONE, config); /*--- Definition of the boundary condition method ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CGeneralAvgGrad_Flow(nDim, nVar_Flow, VISCOUS_GRAD_CORR::NONE, config); } } if (incompressible) { /*--- Incompressible flow, use preconditioning method ---*/ - numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, true, config); + numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, config->GetKind_ViscousGradCorr(), config); for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, VISCOUS_GRAD_CORR::NONE, config); /*--- Definition of the boundary condition method ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, false, config); + numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, VISCOUS_GRAD_CORR::NONE, config); } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -2077,8 +2077,8 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver /*--- Definition of the viscous scheme for each equation and mesh level ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - numerics[iMGlevel][HEAT_SOL][visc_term] = new CAvgGrad_Heat(nDim, config, true); - numerics[iMGlevel][HEAT_SOL][visc_bound_term] = new CAvgGrad_Heat(nDim, config, false); + numerics[iMGlevel][HEAT_SOL][visc_term] = new CAvgGrad_Heat(nDim, config, config->GetKind_ViscousGradCorr()); + numerics[iMGlevel][HEAT_SOL][visc_bound_term] = new CAvgGrad_Heat(nDim, config, VISCOUS_GRAD_CORR::NONE); switch (config->GetKind_ConvNumScheme_Heat()) { diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index ed833df16a2f..8b41136fff09 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -32,7 +32,7 @@ CAvgGrad_Base::CAvgGrad_Base(unsigned short val_nDim, unsigned short val_nVar, unsigned short val_nPrimVar, - bool val_correct_grad, + VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config) : CNumerics(val_nDim, val_nVar, config), nPrimVar(val_nPrimVar), @@ -106,14 +106,52 @@ void CAvgGrad_Base::CorrectGradient(su2double** GradPrimVar, const su2double* val_edge_vector, const su2double val_dist_ij_2, const unsigned short val_nPrimVar) { + + su2double diss[MAXNDIM] = {0.0}; + + const su2double eDotN = 1.0 / GeometryToolbox::DotProduct(nDim,val_edge_vector,UnitNormal); + const su2double alpha = 4.0 / 3.0; + + switch (correct_gradient) { + case VISCOUS_GRAD_CORR::EDGE_NORMAL: + for (int iDim = 0; iDim < nDim; ++iDim) { + diss[iDim] = val_edge_vector[iDim] / val_dist_ij_2; + } + break; + case VISCOUS_GRAD_CORR::FACE_TANGENT: + for (int iDim = 0; iDim < nDim; ++iDim) { + diss[iDim] = UnitNormal[iDim] * eDotN; + } + break; + case VISCOUS_GRAD_CORR::ALPHA_DAMPING: + for (int iDim = 0; iDim < nDim; ++iDim) { + diss[iDim] = UnitNormal[iDim] * alpha * fabs(eDotN); + } + break; + default: + return; + } + + GradientCorrection(GradPrimVar, val_PrimVar_i, val_PrimVar_j, val_edge_vector, + val_dist_ij_2, diss, val_nPrimVar); + +} + +void CAvgGrad_Base::GradientCorrection(su2double** GradPrimVar, + const su2double* val_PrimVar_i, + const su2double* val_PrimVar_j, + const su2double* val_edge_vector, + const su2double val_dist_ij_2, + const su2double* diss, + const unsigned short val_nPrimVar) { for (unsigned short iVar = 0; iVar < val_nPrimVar; iVar++) { Proj_Mean_GradPrimVar_Edge[iVar] = 0.0; for (unsigned short iDim = 0; iDim < nDim; iDim++) { Proj_Mean_GradPrimVar_Edge[iVar] += GradPrimVar[iVar][iDim]*val_edge_vector[iDim]; } for (unsigned short iDim = 0; iDim < nDim; iDim++) { - GradPrimVar[iVar][iDim] -= (Proj_Mean_GradPrimVar_Edge[iVar] - - (val_PrimVar_j[iVar]-val_PrimVar_i[iVar]))*val_edge_vector[iDim] / val_dist_ij_2; + GradPrimVar[iVar][iDim] += diss[iDim] * (val_PrimVar_j[iVar]-val_PrimVar_i[iVar] + - Proj_Mean_GradPrimVar_Edge[iVar]); } } } @@ -363,7 +401,7 @@ void CAvgGrad_Base::GetViscousProjJacs(const su2double *val_Mean_PrimVar, CAvgGrad_Flow::CAvgGrad_Flow(unsigned short val_nDim, unsigned short val_nVar, - bool val_correct_grad, + VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config) : CAvgGrad_Base(val_nDim, val_nVar, val_nDim+3, val_correct_grad, config) { } @@ -429,7 +467,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) /*--- Projection of the mean gradient in the direction of the edge ---*/ - if (correct_gradient && dist_ij_2 != 0.0) { + if (correct_gradient != VISCOUS_GRAD_CORR::NONE && dist_ij_2 != 0.0) { CorrectGradient(Mean_GradPrimVar, PrimVar_i, PrimVar_j, Edge_Vector, dist_ij_2, nDim+1); } @@ -535,7 +573,7 @@ void CAvgGrad_Flow::SetHeatFluxJacobian(const su2double *val_Mean_PrimVar, CAvgGradInc_Flow::CAvgGradInc_Flow(unsigned short val_nDim, unsigned short val_nVar, - bool val_correct_grad, const CConfig* config) + VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config) : CAvgGrad_Base(val_nDim, val_nVar, val_nDim+3, val_correct_grad, config) { energy = config->GetEnergy_Equation(); @@ -600,7 +638,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi /*--- Projection of the mean gradient in the direction of the edge ---*/ - if (correct_gradient && dist_ij_2 != 0.0) { + if (correct_gradient != VISCOUS_GRAD_CORR::NONE && dist_ij_2 != 0.0) { CorrectGradient(Mean_GradPrimVar, PrimVar_i, PrimVar_j, Edge_Vector, dist_ij_2, nVar); } @@ -792,7 +830,7 @@ void CAvgGradInc_Flow::GetViscousIncProjJacs(su2double val_dS, CGeneralAvgGrad_Flow::CGeneralAvgGrad_Flow(unsigned short val_nDim, unsigned short val_nVar, - bool val_correct_grad, + VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config) : CAvgGrad_Base(val_nDim, val_nVar, val_nDim+4, val_correct_grad, config) { } @@ -916,7 +954,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c /*--- Projection of the mean gradient in the direction of the edge ---*/ - if (correct_gradient && dist_ij_2 != 0.0) { + if (correct_gradient != VISCOUS_GRAD_CORR::NONE && dist_ij_2 != 0.0) { CorrectGradient(Mean_GradPrimVar, PrimVar_i, PrimVar_j, Edge_Vector, dist_ij_2, nDim+1); } diff --git a/SU2_CFD/src/numerics/radiation.cpp b/SU2_CFD/src/numerics/radiation.cpp index a52911d3a8a2..d6a067a23567 100644 --- a/SU2_CFD/src/numerics/radiation.cpp +++ b/SU2_CFD/src/numerics/radiation.cpp @@ -89,7 +89,7 @@ void CAvgGradCorrected_P1::ComputeResidual(su2double *val_residual, su2double ** su2double NormalGrad[nVar], CorrectedGrad[nVar]; auto proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, RadVar_Grad_i, - RadVar_Grad_j, true, RadVar_i, RadVar_j, + RadVar_Grad_j, config->GetKind_ViscousGradCorr(), RadVar_i, RadVar_j, NormalGrad, CorrectedGrad); val_residual[0] = GammaP1*CorrectedGrad[0];