From 2a66c5418c9b88def3f52bff4d3a54f823a337ad Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 10 Jun 2024 18:44:12 +0200 Subject: [PATCH 01/44] - Added TKE production limiter constant to config - Added production terms to SST --- Common/include/CConfig.hpp | 7 ++++++- Common/include/option_structure.hpp | 10 ++++++++-- Common/src/CConfig.cpp | 5 +++++ SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 4 ++++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 3 +++ 5 files changed, 26 insertions(+), 3 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 494c9542922..310c0f0f57e 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -856,7 +856,8 @@ class CConfig { nMu_Temperature_Ref, /*!< \brief Number of species reference temperature for Sutherland model. */ nMu_S, /*!< \brief Number of species reference S for Sutherland model. */ nThermal_Conductivity_Constant,/*!< \brief Number of species constant thermal conductivity. */ - nPrandtl_Lam, /*!< \brief Number of species laminar Prandtl number. */ + nPrandtl_Lam, /*!< \brief Number of species + addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); Prandtl number. */ nPrandtl_Turb, /*!< \brief Number of species turbulent Prandtl number. */ nConstant_Lewis_Number; /*!< \brief Number of species Lewis Number. */ su2double Diffusivity_Constant; /*!< \brief Constant mass diffusivity for scalar transport. */ @@ -868,6 +869,8 @@ class CConfig { array MuPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */ array KtPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */ su2double TurbIntensityAndViscRatioFreeStream[2]; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */ + su2double TKE_ProductionLimiterConstant; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */ + bool Change_TKE_ProductionLimiterConstant; su2double Energy_FreeStream, /*!< \brief Free-stream total energy of the fluid. */ ModVel_FreeStream, /*!< \brief Magnitude of the free-stream velocity of the fluid. */ ModVel_FreeStreamND, /*!< \brief Non-dimensional magnitude of the free-stream velocity of the fluid. */ @@ -8813,6 +8816,8 @@ class CConfig { * \brief Set freestream turbonormal for initializing solution. */ const su2double* GetFreeStreamTurboNormal(void) const { return FreeStreamTurboNormal; } + const su2double GetTKE_ProductionLimiterConstant(void) const { return TKE_ProductionLimiterConstant; } + const bool GetChange_TKE_ProductionLimiterConstant(void) const { return Change_TKE_ProductionLimiterConstant; } /*! * \brief Set multizone properties. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 9bbe613c4bc..e355a8364eb 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -990,18 +990,20 @@ enum class SST_OPTIONS { V, /*!< \brief Menter k-w SST model with vorticity production terms. */ KL, /*!< \brief Menter k-w SST model with Kato-Launder production terms. */ UQ, /*!< \brief Menter k-w SST model with uncertainty quantification modifications. */ + FULLPROD, /*!< \brief Menter k-w SST model with full production term. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) MakePair("V1994m", SST_OPTIONS::V1994m) MakePair("V2003m", SST_OPTIONS::V2003m) /// TODO: For now we do not support "unmodified" versions of SST. - //MakePair("V1994", SST_OPTIONS::V1994) - //MakePair("V2003", SST_OPTIONS::V2003) + MakePair("V1994", SST_OPTIONS::V1994) + MakePair("V2003", SST_OPTIONS::V2003) MakePair("SUSTAINING", SST_OPTIONS::SUST) MakePair("VORTICITY", SST_OPTIONS::V) MakePair("KATO-LAUNDER", SST_OPTIONS::KL) MakePair("UQ", SST_OPTIONS::UQ) + MakePair("FULLPROD", SST_OPTIONS::FULLPROD) }; /*! @@ -1013,6 +1015,7 @@ struct SST_ParsedOptions { bool sust = false; /*!< \brief Bool for SST model with sustaining terms. */ bool uq = false; /*!< \brief Bool for using uncertainty quantification. */ bool modified = false; /*!< \brief Bool for modified (m) SST model. */ + bool fullProd = false; /*!< \brief Bool for full production term. */ }; /*! @@ -1030,6 +1033,8 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne return std::find(SST_Options, sst_options_end, option) != sst_options_end; }; + const bool found_fullProd = IsPresent(SST_OPTIONS::FULLPROD); + const bool found_1994 = IsPresent(SST_OPTIONS::V1994); const bool found_2003 = IsPresent(SST_OPTIONS::V2003); const bool found_1994m = IsPresent(SST_OPTIONS::V1994m); @@ -1071,6 +1076,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.sust = sst_sust; SSTParsedOptions.modified = sst_m; SSTParsedOptions.uq = sst_uq; + SSTParsedOptions.fullProd = found_fullProd; return SSTParsedOptions; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 5ff351cfee5..d42634331d3 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1427,6 +1427,8 @@ void CConfig::SetConfig_Options() { addDoubleOption("INITIAL_BCTHRUST", Initial_BCThrust, 4000.0); /* DESCRIPTION: */ addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); + addBoolOption("CHANGE_TKE_PRODUCTION_LIMITER_CONSTANT", Change_TKE_ProductionLimiterConstant, false); + addDoubleOption("TKE_PRODUCTION_LIMITER_CONSTANT", TKE_ProductionLimiterConstant, 10.0); /* DESCRIPTION: Side-slip angle (degrees, only for compressible flows) */ addDoubleOption("SIDESLIP_ANGLE", AoS, 0.0); /*!\brief AOA \n DESCRIPTION: Angle of attack (degrees, only for compressible flows) \ingroup Config*/ @@ -6175,6 +6177,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { break; } cout << "." << endl; + + if (Change_TKE_ProductionLimiterConstant) cout << "Changinf the value of the TKE production limiter constant to " << TKE_ProductionLimiterConstant << endl; + break; } switch (Kind_Trans_Model) { diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index bdf9d0c6fc4..3fa43f57679 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -805,6 +805,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics { const su2double prod_limit = prod_lim_const * beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0]; su2double P = Eddy_Viscosity_i * pow(P_Base, 2); + if (!sstParsedOptions.modified) P -= Eddy_Viscosity_i * diverg*diverg * 2.0/3.0; + if (sstParsedOptions.fullProd) P -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.0; + su2double pk = max(0.0, min(P, prod_limit)); const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; @@ -863,6 +866,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Implicit part ---*/ Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume; + if (sstParsedOptions.fullProd) Jacobian_i[0][0] -= diverg * Volume*2.0/3.0; Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume; Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume; diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 555f8e77155..3a38e5108b8 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -108,6 +108,9 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh constants[9] = 0.44; //gamma_2 constants[10] = 10.0; // production limiter constant } + + if (config->GetChange_TKE_ProductionLimiterConstant()) constants[10] = config->GetTKE_ProductionLimiterConstant(); + /*--- Initialize lower and upper limits---*/ lowerlimit[0] = 1.0e-10; upperlimit[0] = 1.0e10; From 75825906dc03119d520bc7ab58b0aeef60516b0f Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 17 Jun 2024 16:20:50 +0200 Subject: [PATCH 02/44] - Added Production and Destruction terms of SST to output --- SU2_CFD/include/numerics/CNumerics.hpp | 4 +++ .../numerics/turbulent/turb_sources.hpp | 11 +++++++ SU2_CFD/include/variables/CTurbVariable.hpp | 33 +++++++++++++++++++ SU2_CFD/include/variables/CVariable.hpp | 6 ++++ SU2_CFD/src/output/CFlowOutput.cpp | 8 +++++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 7 ++++ SU2_CFD/src/variables/CTurbVariable.cpp | 7 ++++ 7 files changed, 76 insertions(+) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 450112df8f7..5f4545f7b13 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -190,6 +190,8 @@ class CNumerics { bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ + su2double ProdDistr[4]; + public: /*! * \brief Return type used in some "ComputeResidual" overloads to give a @@ -1605,6 +1607,8 @@ class CNumerics { * \return is_bounded_scalar : scalar solver uses bounded scalar convective transport */ inline bool GetBoundedScalar() const { return bounded_scalar;} + + inline su2double GetProdDest(int index) const {return ProdDistr[index];} }; /*! diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 3fa43f57679..f9aeed8f244 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -609,6 +609,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Ambient values for SST-SUST. ---*/ const su2double kAmb, omegaAmb; + su2double Pk, Dk, Pw, Dw; + su2double F1_i, F2_i, CDkw_i; su2double Residual[2]; su2double* Jacobian_i[2]; @@ -846,15 +848,24 @@ class CSourcePieceWise_TurbSST final : public CNumerics { } /*--- Add the production terms to the residuals. ---*/ + Pk = pk; + Pw = pw; Residual[0] += pk * Volume; Residual[1] += pw * Volume; /*--- Add the dissipation terms to the residuals.---*/ + Dk = dk; + Dw = dw; Residual[0] -= dk * Volume; Residual[1] -= dw * Volume; + ProdDistr[0] = Pk; + ProdDistr[1] = Dk; + ProdDistr[2] = Pw; + ProdDistr[3] = Dw; + /*--- Cross diffusion ---*/ Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index eb74e40dbf1..facfac5af65 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -43,6 +43,10 @@ class CTurbVariable : public CScalarVariable { static constexpr size_t MAXNVAR = 2; VectorType turb_index; VectorType intermittency; /*!< \brief Value of the intermittency for the trans. model. */ + VectorType Pk; + VectorType Pw; + VectorType Dk; + VectorType Dw; /*! * \brief Constructor of the class. @@ -100,5 +104,34 @@ class CTurbVariable : public CScalarVariable { */ inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; } + /*! + * \brief Set the intermittency of the transition model. + * \param[in] iPoint - Point index. + * \param[in] val_intermittency - New value of the intermittency. + */ + inline void SetProdDestr(unsigned long iPoint, su2double* val_ProdDestr) final { + Pk(iPoint) = val_ProdDestr[0]; + Dk(iPoint) = val_ProdDestr[1]; + Pw(iPoint) = val_ProdDestr[2]; + Dw(iPoint) = val_ProdDestr[3]; + } + + /*! + * \brief Set the intermittency of the transition model. + * \param[in] iPoint - Point index. + * \param[in] val_intermittency - New value of the intermittency. + */ + inline su2double GetProdTKE(unsigned long iPoint) const final { + return Pk(iPoint); + } + inline su2double GetDestrTKE(unsigned long iPoint) const final { + return Dk(iPoint); + } + inline su2double GetProdW(unsigned long iPoint) const final { + return Pw(iPoint); + } + inline su2double GetDestrW(unsigned long iPoint) const final { + return Dw(iPoint); + } }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 358824ae8d7..2d340ee0c39 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2340,4 +2340,10 @@ class CVariable { inline virtual const su2double *GetScalarSources(unsigned long iPoint) const { return nullptr; } inline virtual const su2double *GetScalarLookups(unsigned long iPoint) const { return nullptr; } + + inline virtual void SetProdDestr(unsigned long iPoint, su2double* val_ProdDestr) { } + inline virtual su2double GetProdTKE(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetDestrTKE(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetProdW(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetDestrW(unsigned long iPoint) const { return 0.0; } }; diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 5ffee0d23f5..b0150173ac8 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1252,6 +1252,10 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ case TURB_FAMILY::KW: AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy"); AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); + AddVolumeOutput("PROD_TKE", "Prod_TKE", "SOLUTION", "Production of turbulent kinetic energy"); + AddVolumeOutput("DESTR_TKE", "Destr_TKE", "SOLUTION", "Destruction of turbulent kinetic energy"); + AddVolumeOutput("PROD_W", "Prod_W", "SOLUTION", "Production of rate of dissipation"); + AddVolumeOutput("DESTR_W", "Destr_W", "SOLUTION", "Destruction of rate of dissipation"); break; case TURB_FAMILY::NONE: @@ -1528,6 +1532,10 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con case TURB_FAMILY::KW: SetVolumeOutputValue("TKE", iPoint, Node_Turb->GetSolution(iPoint, 0)); SetVolumeOutputValue("DISSIPATION", iPoint, Node_Turb->GetSolution(iPoint, 1)); + SetVolumeOutputValue("PROD_TKE", iPoint, Node_Turb->GetProdTKE(iPoint)); + SetVolumeOutputValue("DESTR_TKE", iPoint, Node_Turb->GetDestrTKE(iPoint)); + SetVolumeOutputValue("PROD_W", iPoint, Node_Turb->GetProdW(iPoint)); + SetVolumeOutputValue("DESTR_W", iPoint, Node_Turb->GetDestrW(iPoint)); SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); if (limiter) { diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 3a38e5108b8..cda1a677904 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -371,6 +371,13 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta nodes->SetIntermittency(iPoint, numerics->GetIntermittencyEff()); } + su2double ProdDestr[4]; + ProdDestr[0] = numerics->GetProdDest(0); + ProdDestr[1] = numerics->GetProdDest(1); + ProdDestr[2] = numerics->GetProdDest(2); + ProdDestr[3] = numerics->GetProdDest(3); + nodes->SetProdDestr(iPoint, ProdDestr); + /*--- Subtract residual and the Jacobian ---*/ LinSysRes.SubtractBlock(iPoint, residual); diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp index 7fcc357c0ee..773bb78dd1d 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -35,4 +35,11 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned turb_index.resize(nPoint) = su2double(1.0); intermittency.resize(nPoint) = su2double(1.0); + if (TurbModelFamily(config->GetKind_Turb_Model()) == TURB_FAMILY::KW) { + Pk.resize(nPoint) = su2double(1.0); + Pw.resize(nPoint) = su2double(1.0); + Dk.resize(nPoint) = su2double(1.0); + Dw.resize(nPoint) = su2double(1.0); + } + } From edce8ee37533cf6b86dc0180dadd0c85fb749a5c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 19 Jul 2024 17:12:13 +0200 Subject: [PATCH 03/44] - Added user-defined lower limits of TKE and W - Added user-defined production limiter constant for TKE --- Common/include/CConfig.hpp | 7 +++++++ Common/include/option_structure.hpp | 10 ++++++++++ Common/src/CConfig.cpp | 4 ++++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 7 +++++++ 4 files changed, 28 insertions(+) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 310c0f0f57e..042a866ed56 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1168,6 +1168,9 @@ class CConfig { nHistoryOutput, nVolumeOutput; /*!< \brief Number of variables printed to the history file. */ bool Multizone_Residual; /*!< \brief Determines if memory should be allocated for the multizone residual. */ SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */ + su2double lowerLimitTKE, + lowerLimitDissipation; + su2double prodLimConst; SA_ParsedOptions saParsedOptions; /*!< \brief Additional parameters for the SA turbulence model. */ LM_ParsedOptions lmParsedOptions; /*!< \brief Additional parameters for the LM transition model. */ su2double uq_delta_b; /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */ @@ -9852,6 +9855,10 @@ class CConfig { */ SST_ParsedOptions GetSSTParsedOptions() const { return sstParsedOptions; } + su2double GetLowerLimitTKE() const { return lowerLimitTKE; } + su2double GetLowerLimitDissipation() const { return lowerLimitDissipation; } + su2double GetProdLimConst() const { return prodLimConst; } + /*! * \brief Get parsed SA option data structure. * \return SA option data structure. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index e355a8364eb..4b5e0b0c071 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -991,6 +991,8 @@ enum class SST_OPTIONS { KL, /*!< \brief Menter k-w SST model with Kato-Launder production terms. */ UQ, /*!< \brief Menter k-w SST model with uncertainty quantification modifications. */ FULLPROD, /*!< \brief Menter k-w SST model with full production term. */ + LLT, /*!< \brief Menter k-w SST model with Lower Limits Treatments. */ + PRODLIM, /*!< \brief Menter k-w SST model with user-defined production limiter constant. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) @@ -1004,6 +1006,8 @@ static const MapType SST_Options_Map = { MakePair("KATO-LAUNDER", SST_OPTIONS::KL) MakePair("UQ", SST_OPTIONS::UQ) MakePair("FULLPROD", SST_OPTIONS::FULLPROD) + MakePair("LLT", SST_OPTIONS::LLT) + MakePair("PRODLIM", SST_OPTIONS::PRODLIM) }; /*! @@ -1016,6 +1020,8 @@ struct SST_ParsedOptions { bool uq = false; /*!< \brief Bool for using uncertainty quantification. */ bool modified = false; /*!< \brief Bool for modified (m) SST model. */ bool fullProd = false; /*!< \brief Bool for full production term. */ + bool llt = false; /*!< \brief Bool for Lower Limits Treatment. */ + bool prodLim = false; /*!< \brief Bool for user-defined production limiter constant. */ }; /*! @@ -1034,6 +1040,8 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne }; const bool found_fullProd = IsPresent(SST_OPTIONS::FULLPROD); + const bool found_llt = IsPresent(SST_OPTIONS::LLT); + const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM); const bool found_1994 = IsPresent(SST_OPTIONS::V1994); const bool found_2003 = IsPresent(SST_OPTIONS::V2003); @@ -1077,6 +1085,8 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.modified = sst_m; SSTParsedOptions.uq = sst_uq; SSTParsedOptions.fullProd = found_fullProd; + SSTParsedOptions.llt = found_llt; + SSTParsedOptions.prodLim = found_prodLim; return SSTParsedOptions; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index d42634331d3..8297075bbe3 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1118,6 +1118,10 @@ void CConfig::SetConfig_Options() { /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SA turbulence model options/corrections. \n Options: see \link SA_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumListOption("SA_OPTIONS", nSA_Options, SA_Options, SA_Options_Map); + addDoubleOption("LOWER_LIMIT_TKE", lowerLimitTKE, 1e-20); + addDoubleOption("LOWER_LIMIT_DISSIPATION", lowerLimitDissipation, 1e-6); + addDoubleOption("PROD_LIM_CONST", prodLimConst, 20.0); + /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumOption("KIND_TRANS_MODEL", Kind_Trans_Model, Trans_Model_Map, TURB_TRANS_MODEL::NONE); /*!\brief SST_OPTIONS \n DESCRIPTION: Specify LM transition model options/correlations. \n Options: see \link LM_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index cda1a677904..9a1f6647e66 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -102,11 +102,13 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh constants[8] = constants[4]/constants[6] - constants[2]*0.41*0.41/sqrt(constants[6]); //alfa_1 constants[9] = constants[5]/constants[6] - constants[3]*0.41*0.41/sqrt(constants[6]); //alfa_2 constants[10] = 20.0; // production limiter constant + if (sstParsedOptions.prodLim) constants[10] = config->GetProdLimConst(); } else { /* SST-V2003 */ constants[8] = 5.0 / 9.0; //gamma_1 constants[9] = 0.44; //gamma_2 constants[10] = 10.0; // production limiter constant + if (sstParsedOptions.prodLim) constants[10] = config->GetProdLimConst(); } if (config->GetChange_TKE_ProductionLimiterConstant()) constants[10] = config->GetTKE_ProductionLimiterConstant(); @@ -135,6 +137,11 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh Solution_Inf[0] = kine_Inf; Solution_Inf[1] = omega_Inf; + if (sstParsedOptions.llt) { + lowerlimit[0] = kine_Inf * config->GetLowerLimitTKE(); + lowerlimit[1] = omega_Inf * config->GetLowerLimitDissipation(); + } + /*--- Eddy viscosity, initialized without stress limiter at the infinity ---*/ muT_Inf = rhoInf*kine_Inf/omega_Inf; From 41c61cac187d6fae0490b9e6f80a8bdac962d0a7 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 22 Jul 2024 09:50:27 +0200 Subject: [PATCH 04/44] - Changed the Cross-Diffusion term in F1 computations --- SU2_CFD/src/variables/CTurbSSTVariable.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index 3b88c4b8f3f..eae2c3ef0b2 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -69,7 +69,10 @@ void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_visco for (unsigned long iDim = 0; iDim < nDim; iDim++) CDkw(iPoint) += Gradient(iPoint,0,iDim)*Gradient(iPoint,1,iDim); CDkw(iPoint) *= 2.0*val_density*sigma_om2/Solution(iPoint,1); - CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -prod_lim_const)); + // CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -prod_lim_const)); + su2double exponent = 20.0; + if (sstParsedOptions.version == SST_OPTIONS::V2003) exponent = 10.0; + CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -exponent)); /*--- F1 ---*/ From 7b11baa0b2b766458f67c6f23c137154caccc5ff Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 22 Jul 2024 09:56:30 +0200 Subject: [PATCH 05/44] - Removed obsolete changes to TKE prod limiter --- Common/include/CConfig.hpp | 7 ++----- Common/src/CConfig.cpp | 2 -- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 1 - 3 files changed, 2 insertions(+), 8 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 042a866ed56..e3ce7a2dcab 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -869,8 +869,7 @@ class CConfig { array MuPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */ array KtPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */ su2double TurbIntensityAndViscRatioFreeStream[2]; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */ - su2double TKE_ProductionLimiterConstant; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */ - bool Change_TKE_ProductionLimiterConstant; + su2double Energy_FreeStream, /*!< \brief Free-stream total energy of the fluid. */ ModVel_FreeStream, /*!< \brief Magnitude of the free-stream velocity of the fluid. */ ModVel_FreeStreamND, /*!< \brief Non-dimensional magnitude of the free-stream velocity of the fluid. */ @@ -8819,8 +8818,6 @@ class CConfig { * \brief Set freestream turbonormal for initializing solution. */ const su2double* GetFreeStreamTurboNormal(void) const { return FreeStreamTurboNormal; } - const su2double GetTKE_ProductionLimiterConstant(void) const { return TKE_ProductionLimiterConstant; } - const bool GetChange_TKE_ProductionLimiterConstant(void) const { return Change_TKE_ProductionLimiterConstant; } /*! * \brief Set multizone properties. @@ -9858,7 +9855,7 @@ class CConfig { su2double GetLowerLimitTKE() const { return lowerLimitTKE; } su2double GetLowerLimitDissipation() const { return lowerLimitDissipation; } su2double GetProdLimConst() const { return prodLimConst; } - + /*! * \brief Get parsed SA option data structure. * \return SA option data structure. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 8297075bbe3..943acf94759 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1431,8 +1431,6 @@ void CConfig::SetConfig_Options() { addDoubleOption("INITIAL_BCTHRUST", Initial_BCThrust, 4000.0); /* DESCRIPTION: */ addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); - addBoolOption("CHANGE_TKE_PRODUCTION_LIMITER_CONSTANT", Change_TKE_ProductionLimiterConstant, false); - addDoubleOption("TKE_PRODUCTION_LIMITER_CONSTANT", TKE_ProductionLimiterConstant, 10.0); /* DESCRIPTION: Side-slip angle (degrees, only for compressible flows) */ addDoubleOption("SIDESLIP_ANGLE", AoS, 0.0); /*!\brief AOA \n DESCRIPTION: Angle of attack (degrees, only for compressible flows) \ingroup Config*/ diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 9a1f6647e66..b0b5cabcd7f 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -111,7 +111,6 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh if (sstParsedOptions.prodLim) constants[10] = config->GetProdLimConst(); } - if (config->GetChange_TKE_ProductionLimiterConstant()) constants[10] = config->GetTKE_ProductionLimiterConstant(); /*--- Initialize lower and upper limits---*/ lowerlimit[0] = 1.0e-10; From 9482f513ae58e7ce0a3bf3a67c6783de75232d98 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 22 Jul 2024 09:57:48 +0200 Subject: [PATCH 06/44] - Modified the config output --- Common/src/CConfig.cpp | 2 +- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 943acf94759..50d70692889 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6180,7 +6180,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { } cout << "." << endl; - if (Change_TKE_ProductionLimiterConstant) cout << "Changinf the value of the TKE production limiter constant to " << TKE_ProductionLimiterConstant << endl; + if (sstParsedOptions.prodLim) cout << "Changing the value of the TKE production limiter constant to " << prodLimConst << endl; break; } diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index b0b5cabcd7f..7d4b5d4e430 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -110,7 +110,6 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh constants[10] = 10.0; // production limiter constant if (sstParsedOptions.prodLim) constants[10] = config->GetProdLimConst(); } - /*--- Initialize lower and upper limits---*/ lowerlimit[0] = 1.0e-10; From 38a02d63306d83338b5d5965e20d2d578aad1a45 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 22 Jul 2024 10:00:20 +0200 Subject: [PATCH 07/44] - Added lower limit changes for SST in config output --- Common/src/CConfig.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 50d70692889..96bdf54e6f3 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6181,6 +6181,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "." << endl; if (sstParsedOptions.prodLim) cout << "Changing the value of the TKE production limiter constant to " << prodLimConst << endl; + if (sstParsedOptions.llt) cout << "Changing the value of the lower limits of TKE and Omega " << endl; break; } From c3e0d4b8cfd240522db6b284ee9dce31c5daf9b9 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 22 Jul 2024 10:04:09 +0200 Subject: [PATCH 08/44] - Added wall distance output for mesh adaptation --- SU2_CFD/src/output/CFlowOutput.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index b0150173ac8..01aabe0fee2 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1483,6 +1483,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) { AddVolumeOutput("VORTICITY", "Vorticity", "VORTEX_IDENTIFICATION", "Value of the vorticity"); } AddVolumeOutput("Q_CRITERION", "Q_Criterion", "VORTEX_IDENTIFICATION", "Value of the Q-Criterion"); + AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "MISC", "Wall distance value"); } // Timestep info @@ -1516,6 +1517,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]); } SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(Node_Flow->GetVelocityGradient(iPoint))); + SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint)); } const bool limiter = (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE); From 0f7558b40c6c94781ebf42c23f1b098cef132e5d Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 22 Jul 2024 11:12:21 +0200 Subject: [PATCH 09/44] - Added strain magnitude as output --- SU2_CFD/src/output/CFlowOutput.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 01aabe0fee2..cf614835a53 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1484,6 +1484,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) { } AddVolumeOutput("Q_CRITERION", "Q_Criterion", "VORTEX_IDENTIFICATION", "Value of the Q-Criterion"); AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "MISC", "Wall distance value"); + AddVolumeOutput("STRAIN_MAG", "Strain_Magnitude", "MISC", "Strain magnitude value"); } // Timestep info @@ -1518,6 +1519,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con } SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(Node_Flow->GetVelocityGradient(iPoint))); SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint)); + SetVolumeOutputValue("STRAIN_MAG", iPoint, Node_Turb->GetStrainMag(iPoint)); } const bool limiter = (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE); From 3be370d8e44337c4e6b1d64dfb6e4e59558b7553 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 22 Jul 2024 11:56:20 +0200 Subject: [PATCH 10/44] - added production limiter flag to output --- SU2_CFD/include/numerics/CNumerics.hpp | 2 +- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 3 +++ SU2_CFD/include/variables/CTurbVariable.hpp | 5 +++++ SU2_CFD/include/variables/CVariable.hpp | 1 + SU2_CFD/src/output/CFlowOutput.cpp | 2 ++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 3 ++- SU2_CFD/src/variables/CTurbVariable.cpp | 1 + 7 files changed, 15 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 5f4545f7b13..151bcdc2404 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -190,7 +190,7 @@ class CNumerics { bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ - su2double ProdDistr[4]; + su2double ProdDistr[5]; public: /*! diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index f9aeed8f244..482f0d4b089 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -810,6 +810,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics { if (!sstParsedOptions.modified) P -= Eddy_Viscosity_i * diverg*diverg * 2.0/3.0; if (sstParsedOptions.fullProd) P -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.0; + su2double PLim = 0.0; + if ( P > prod_limit ) PLim = 1.0; su2double pk = max(0.0, min(P, prod_limit)); const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; @@ -865,6 +867,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { ProdDistr[1] = Dk; ProdDistr[2] = Pw; ProdDistr[3] = Dw; + ProdDistr[4] = PLim; /*--- Cross diffusion ---*/ diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index facfac5af65..a95ca9f0991 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -47,6 +47,7 @@ class CTurbVariable : public CScalarVariable { VectorType Pw; VectorType Dk; VectorType Dw; + VectorType PkLim; /*! * \brief Constructor of the class. @@ -114,6 +115,7 @@ class CTurbVariable : public CScalarVariable { Dk(iPoint) = val_ProdDestr[1]; Pw(iPoint) = val_ProdDestr[2]; Dw(iPoint) = val_ProdDestr[3]; + PkLim(iPoint) = val_ProdDestr[4]; } /*! @@ -133,5 +135,8 @@ class CTurbVariable : public CScalarVariable { inline su2double GetDestrW(unsigned long iPoint) const final { return Dw(iPoint); } + inline su2double GetPkLim(unsigned long iPoint) const final { + return PkLim(iPoint); + } }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 2d340ee0c39..0f8360c18c9 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2346,4 +2346,5 @@ class CVariable { inline virtual su2double GetDestrTKE(unsigned long iPoint) const { return 0.0; } inline virtual su2double GetProdW(unsigned long iPoint) const { return 0.0; } inline virtual su2double GetDestrW(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetPkLim(unsigned long iPoint) const { return 0.0; } }; diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index cf614835a53..f5607c62965 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1254,6 +1254,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); AddVolumeOutput("PROD_TKE", "Prod_TKE", "SOLUTION", "Production of turbulent kinetic energy"); AddVolumeOutput("DESTR_TKE", "Destr_TKE", "SOLUTION", "Destruction of turbulent kinetic energy"); + AddVolumeOutput("PROD_TKE_LIM", "Prod_TKE_Lim", "SOLUTION", "Check if production limiter has been used for TKE"); AddVolumeOutput("PROD_W", "Prod_W", "SOLUTION", "Production of rate of dissipation"); AddVolumeOutput("DESTR_W", "Destr_W", "SOLUTION", "Destruction of rate of dissipation"); break; @@ -1538,6 +1539,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("DISSIPATION", iPoint, Node_Turb->GetSolution(iPoint, 1)); SetVolumeOutputValue("PROD_TKE", iPoint, Node_Turb->GetProdTKE(iPoint)); SetVolumeOutputValue("DESTR_TKE", iPoint, Node_Turb->GetDestrTKE(iPoint)); + SetVolumeOutputValue("PROD_TKE_LIM", iPoint, Node_Turb->GetPkLim(iPoint)); SetVolumeOutputValue("PROD_W", iPoint, Node_Turb->GetProdW(iPoint)); SetVolumeOutputValue("DESTR_W", iPoint, Node_Turb->GetDestrW(iPoint)); SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 7d4b5d4e430..74c3a628aa2 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -376,11 +376,12 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta nodes->SetIntermittency(iPoint, numerics->GetIntermittencyEff()); } - su2double ProdDestr[4]; + su2double ProdDestr[5]; ProdDestr[0] = numerics->GetProdDest(0); ProdDestr[1] = numerics->GetProdDest(1); ProdDestr[2] = numerics->GetProdDest(2); ProdDestr[3] = numerics->GetProdDest(3); + ProdDestr[4] = numerics->GetProdDest(4); nodes->SetProdDestr(iPoint, ProdDestr); /*--- Subtract residual and the Jacobian ---*/ diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp index 773bb78dd1d..00c9d7b2fa0 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -40,6 +40,7 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned Pw.resize(nPoint) = su2double(1.0); Dk.resize(nPoint) = su2double(1.0); Dw.resize(nPoint) = su2double(1.0); + PkLim.resize(nPoint) = su2double(1.0); } } From 125fd6df50ae619eddb7237ad51524e31928b06a Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 24 Jul 2024 15:56:50 +0200 Subject: [PATCH 11/44] - Changed the Cross-Diffusion term in the residual of Omega --- SU2_CFD/src/variables/CTurbSSTVariable.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index eae2c3ef0b2..de200423ca6 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -72,14 +72,14 @@ void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_visco // CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -prod_lim_const)); su2double exponent = 20.0; if (sstParsedOptions.version == SST_OPTIONS::V2003) exponent = 10.0; - CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -exponent)); + const su2double CDkw_Clipped = max(CDkw(iPoint), pow(10.0, -exponent)); /*--- F1 ---*/ arg2A = sqrt(Solution(iPoint,0))/(beta_star*Solution(iPoint,1)*val_dist+EPS*EPS); arg2B = 500.0*val_viscosity / (val_density*val_dist*val_dist*Solution(iPoint,1)+EPS*EPS); arg2 = max(arg2A, arg2B); - arg1 = min(arg2, 4.0*val_density*sigma_om2*Solution(iPoint,0) / (CDkw(iPoint)*val_dist*val_dist+EPS*EPS)); + arg1 = min(arg2, 4.0*val_density*sigma_om2*Solution(iPoint,0) / (CDkw_Clipped*val_dist*val_dist+EPS*EPS)); F1(iPoint) = tanh(pow(arg1, 4.0)); /*--- F2 ---*/ @@ -94,7 +94,7 @@ void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_visco F1(iPoint) = max(F1(iPoint), F3); } - AD::SetPreaccOut(F1(iPoint)); AD::SetPreaccOut(F2(iPoint)); AD::SetPreaccOut(CDkw(iPoint)); + AD::SetPreaccOut(F1(iPoint)); AD::SetPreaccOut(F2(iPoint)); AD::SetPreaccOut(CDkw_Clipped); AD::EndPreacc(); } From 91cc3b8abe9eb5a03543e4c5b1ff75d2642842e3 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 24 Jul 2024 16:29:48 +0200 Subject: [PATCH 12/44] - Added clipping of the cross-diffusion term in SST residual --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 482f0d4b089..80e85f33a3d 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -870,8 +870,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics { ProdDistr[4] = PLim; /*--- Cross diffusion ---*/ - - Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; + su2double CrossDiffusionTerm = (1.0 - F1_i) * CDkw_i; + if (CrossDiffusionTerm < 0.0) CrossDiffusionTerm = -min(Pw, abs(CrossDiffusionTerm)); + Residual[1] += CrossDiffusionTerm * Volume; /*--- Contribution due to 2D axisymmetric formulation ---*/ From 7b78756f530892af36350430d1852b76a20a2c7d Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 24 Jul 2024 17:07:33 +0200 Subject: [PATCH 13/44] - changed cross diffusion term clipping in Omega residual --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 80e85f33a3d..fe12364ff52 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -870,9 +870,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics { ProdDistr[4] = PLim; /*--- Cross diffusion ---*/ - su2double CrossDiffusionTerm = (1.0 - F1_i) * CDkw_i; - if (CrossDiffusionTerm < 0.0) CrossDiffusionTerm = -min(Pw, abs(CrossDiffusionTerm)); - Residual[1] += CrossDiffusionTerm * Volume; + const su2double CrossDiffusionTerm = CDkw_i < 0.0 ? -min(Pw, abs(CrossDiffusionTerm)) : CrossDiffusionTerm; + Residual[1] += (1.0 - F1_i) * CrossDiffusionTerm * Volume; /*--- Contribution due to 2D axisymmetric formulation ---*/ From 076d471fb0d093873b6e2aa672ab00fe8023cc08 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 24 Jul 2024 17:22:07 +0200 Subject: [PATCH 14/44] - Added CDkw and F1 as output --- SU2_CFD/src/output/CFlowOutput.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index f5607c62965..a55134a5066 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1257,6 +1257,8 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ AddVolumeOutput("PROD_TKE_LIM", "Prod_TKE_Lim", "SOLUTION", "Check if production limiter has been used for TKE"); AddVolumeOutput("PROD_W", "Prod_W", "SOLUTION", "Production of rate of dissipation"); AddVolumeOutput("DESTR_W", "Destr_W", "SOLUTION", "Destruction of rate of dissipation"); + AddVolumeOutput("CDkw", "CDkw", "SOLUTION", "Cross-Diffusion term"); + AddVolumeOutput("F1", "F1", "SOLUTION", "F1 blending function"); break; case TURB_FAMILY::NONE: @@ -1542,6 +1544,8 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("PROD_TKE_LIM", iPoint, Node_Turb->GetPkLim(iPoint)); SetVolumeOutputValue("PROD_W", iPoint, Node_Turb->GetProdW(iPoint)); SetVolumeOutputValue("DESTR_W", iPoint, Node_Turb->GetDestrW(iPoint)); + SetVolumeOutputValue("CDkw", iPoint, Node_Turb->GetCrossDiff(iPoint)); + SetVolumeOutputValue("F1", iPoint, Node_Turb->GetF1blending(iPoint)); SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); if (limiter) { From 2e73f775073f1c66a88d28956818fe90d1b4871b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 25 Jul 2024 15:25:29 +0200 Subject: [PATCH 15/44] - Fixed Supersonic inlet BC inclusion of TKE --- SU2_CFD/src/solvers/CEulerSolver.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 356c759564c..6ab73156fd9 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -7341,6 +7341,8 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); /*--- Supersonic inlet flow: there are no outgoing characteristics, so all flow variables can be imposed at the inlet. @@ -7361,14 +7363,15 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con /*--- Compute the energy from the specified state ---*/ const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity); - su2double Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); - + const su2double Energy_woTKE = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; + su2double Energy = Energy_woTKE; + /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + if (tkeNeeded) Energy = Energy_woTKE + turbNodes->GetSolution(iPoint,0); if (!geometry->nodes->GetDomain(iPoint)) continue; From c662daf3d82af32fa19f9a610ee5d7866fbe0323 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 25 Jul 2024 15:41:20 +0200 Subject: [PATCH 16/44] - Corrected Riemann BC for TKE --- SU2_CFD/src/solvers/CEulerSolver.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 6ab73156fd9..b0944bbb036 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -5011,6 +5011,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT), gravity = config->GetGravityForce(), tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); su2double **P_Tensor = new su2double*[nVar], **invP_Tensor = new su2double*[nVar]; @@ -5110,7 +5112,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); + // if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0); break; case STATIC_SUPERSONIC_INFLOW_PT: @@ -7365,7 +7368,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity); const su2double Energy_woTKE = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; su2double Energy = Energy_woTKE; - + /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) From 41e1c16ff9029c4b2244c3f303eb830f77c7a1ed Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 25 Jul 2024 19:41:02 +0200 Subject: [PATCH 17/44] - Fixed more BCs - Updated boundary conditions as in TMR page --- Common/include/CConfig.hpp | 2 ++ Common/include/option_structure.hpp | 5 +++++ Common/src/CConfig.cpp | 1 + SU2_CFD/src/solvers/CEulerSolver.cpp | 9 ++++++++- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 16 ++++++++++++++-- 5 files changed, 30 insertions(+), 3 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index e3ce7a2dcab..aecc1940207 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1170,6 +1170,7 @@ class CConfig { su2double lowerLimitTKE, lowerLimitDissipation; su2double prodLimConst; + su2double LDomain; SA_ParsedOptions saParsedOptions; /*!< \brief Additional parameters for the SA turbulence model. */ LM_ParsedOptions lmParsedOptions; /*!< \brief Additional parameters for the LM transition model. */ su2double uq_delta_b; /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */ @@ -9855,6 +9856,7 @@ class CConfig { su2double GetLowerLimitTKE() const { return lowerLimitTKE; } su2double GetLowerLimitDissipation() const { return lowerLimitDissipation; } su2double GetProdLimConst() const { return prodLimConst; } + su2double GetLDomain() const { return LDomain; } /*! * \brief Get parsed SA option data structure. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 4b5e0b0c071..bec6739be30 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -993,6 +993,7 @@ enum class SST_OPTIONS { FULLPROD, /*!< \brief Menter k-w SST model with full production term. */ LLT, /*!< \brief Menter k-w SST model with Lower Limits Treatments. */ PRODLIM, /*!< \brief Menter k-w SST model with user-defined production limiter constant. */ + NEWBC, /*!< \brief Menter k-w SST model with new boundary conditions. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) @@ -1008,6 +1009,7 @@ static const MapType SST_Options_Map = { MakePair("FULLPROD", SST_OPTIONS::FULLPROD) MakePair("LLT", SST_OPTIONS::LLT) MakePair("PRODLIM", SST_OPTIONS::PRODLIM) + MakePair("NEWBC", SST_OPTIONS::NEWBC) }; /*! @@ -1022,6 +1024,7 @@ struct SST_ParsedOptions { bool fullProd = false; /*!< \brief Bool for full production term. */ bool llt = false; /*!< \brief Bool for Lower Limits Treatment. */ bool prodLim = false; /*!< \brief Bool for user-defined production limiter constant. */ + bool newBC = false; /*!< \brief Bool for new boundary conditions. */ }; /*! @@ -1042,6 +1045,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne const bool found_fullProd = IsPresent(SST_OPTIONS::FULLPROD); const bool found_llt = IsPresent(SST_OPTIONS::LLT); const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM); + const bool found_newBC = IsPresent(SST_OPTIONS::NEWBC); const bool found_1994 = IsPresent(SST_OPTIONS::V1994); const bool found_2003 = IsPresent(SST_OPTIONS::V2003); @@ -1087,6 +1091,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.fullProd = found_fullProd; SSTParsedOptions.llt = found_llt; SSTParsedOptions.prodLim = found_prodLim; + SSTParsedOptions.newBC = found_newBC; return SSTParsedOptions; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 96bdf54e6f3..0c27bc9a49d 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1121,6 +1121,7 @@ void CConfig::SetConfig_Options() { addDoubleOption("LOWER_LIMIT_TKE", lowerLimitTKE, 1e-20); addDoubleOption("LOWER_LIMIT_DISSIPATION", lowerLimitDissipation, 1e-6); addDoubleOption("PROD_LIM_CONST", prodLimConst, 20.0); + addDoubleOption("L_DOMAIN", LDomain, 1.0); /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumOption("KIND_TRANS_MODEL", Kind_Trans_Model, Trans_Model_Map, TURB_TRANS_MODEL::NONE); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index b0944bbb036..ff0227e61b5 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -4768,6 +4768,8 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; bool viscous = config->GetViscous(); bool tkeNeeded = config->GetKind_Turb_Model() == TURB_MODEL::SST; + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); auto *Normal = new su2double[nDim]; @@ -4917,7 +4919,8 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, } Pressure = Density*SoundSpeed*SoundSpeed/Gamma; Energy = Pressure/(Gamma_Minus_One*Density) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + // if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Store new primitive state for computing the flux. ---*/ @@ -7179,6 +7182,9 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, bool gravity = (config->GetGravityForce()); bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); + auto *Normal = new su2double[nDim]; /*--- Loop over all the vertices on this boundary marker ---*/ @@ -7261,6 +7267,7 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, } Energy = P_Exit/(Density*Gamma_Minus_One) + 0.5*Velocity2; if (tkeNeeded) Energy += GetTke_Inf(); + // if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Conservative variables, using the derived quantities ---*/ V_outlet[0] = Pressure / ( Gas_Constant * Density); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 74c3a628aa2..3fb642f1dfc 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -132,6 +132,11 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh su2double kine_Inf = 3.0/2.0*(VelMag2*Intensity*Intensity); su2double omega_Inf = rhoInf*kine_Inf/(muLamInf*viscRatio); + if (sstParsedOptions.newBC) { + su2double omega_Inf = 10 * sqrt(VelMag2) / config->GetLDomain(); + su2double kine_Inf = omega_Inf*(muLamInf*viscRatio)/rhoInf; + } + Solution_Inf[0] = kine_Inf; Solution_Inf[1] = omega_Inf; @@ -632,8 +637,15 @@ void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C const su2double viscRatio = Turb_Properties[1]; const su2double VelMag2 = GeometryToolbox::SquaredNorm(nDim, Velocity_Inlet); - Inlet_Vars[0] = 3.0 / 2.0 * (VelMag2 * pow(Intensity, 2)); - Inlet_Vars[1] = Density_Inlet * Inlet_Vars[0] / (Laminar_Viscosity_Inlet * viscRatio); + if (sstParsedOptions.newBC) { + su2double LDomain = 0.53; + Inlet_Vars[1] = 10 * sqrt(VelMag2) / LDomain; + Inlet_Vars[0] = Inlet_Vars[1]*(Laminar_Viscosity_Inlet*viscRatio)/Density_Inlet; + } else { + Inlet_Vars[0] = 3.0 / 2.0 * (VelMag2 * pow(Intensity, 2)); + Inlet_Vars[1] = Density_Inlet * Inlet_Vars[0] / (Laminar_Viscosity_Inlet * viscRatio); + } + } /*--- Set the turbulent variable states. Use free-stream SST From e9dfc27eef51714f6c0e73babf7d4b525e57a4cc Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 26 Jul 2024 12:48:18 +0200 Subject: [PATCH 18/44] - Fixed bug with Cross diffusion in W residual - Given option for cross diffusion limiting in W residual --- Common/include/CConfig.hpp | 3 +-- Common/include/option_structure.hpp | 5 +++++ SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 10 +++++++++- 3 files changed, 15 insertions(+), 3 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index aecc1940207..8be15998e16 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -856,8 +856,7 @@ class CConfig { nMu_Temperature_Ref, /*!< \brief Number of species reference temperature for Sutherland model. */ nMu_S, /*!< \brief Number of species reference S for Sutherland model. */ nThermal_Conductivity_Constant,/*!< \brief Number of species constant thermal conductivity. */ - nPrandtl_Lam, /*!< \brief Number of species - addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); Prandtl number. */ + nPrandtl_Lam, /*!< \brief Prandtl number. */ nPrandtl_Turb, /*!< \brief Number of species turbulent Prandtl number. */ nConstant_Lewis_Number; /*!< \brief Number of species Lewis Number. */ su2double Diffusivity_Constant; /*!< \brief Constant mass diffusivity for scalar transport. */ diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bec6739be30..3ba4fcf3621 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -994,6 +994,7 @@ enum class SST_OPTIONS { LLT, /*!< \brief Menter k-w SST model with Lower Limits Treatments. */ PRODLIM, /*!< \brief Menter k-w SST model with user-defined production limiter constant. */ NEWBC, /*!< \brief Menter k-w SST model with new boundary conditions. */ + LCD, /*!< \brief Menter k-w SST model with lower limit on the cross-diffusion term. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) @@ -1010,6 +1011,7 @@ static const MapType SST_Options_Map = { MakePair("LLT", SST_OPTIONS::LLT) MakePair("PRODLIM", SST_OPTIONS::PRODLIM) MakePair("NEWBC", SST_OPTIONS::NEWBC) + MakePair("LCD", SST_OPTIONS::LCD) }; /*! @@ -1025,6 +1027,7 @@ struct SST_ParsedOptions { bool llt = false; /*!< \brief Bool for Lower Limits Treatment. */ bool prodLim = false; /*!< \brief Bool for user-defined production limiter constant. */ bool newBC = false; /*!< \brief Bool for new boundary conditions. */ + bool lcd = true; /*!< \brief Bool lower limit on the cross-diffusion term. */ }; /*! @@ -1046,6 +1049,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne const bool found_llt = IsPresent(SST_OPTIONS::LLT); const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM); const bool found_newBC = IsPresent(SST_OPTIONS::NEWBC); + const bool found_lcd = IsPresent(SST_OPTIONS::LCD); const bool found_1994 = IsPresent(SST_OPTIONS::V1994); const bool found_2003 = IsPresent(SST_OPTIONS::V2003); @@ -1092,6 +1096,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.llt = found_llt; SSTParsedOptions.prodLim = found_prodLim; SSTParsedOptions.newBC = found_newBC; + SSTParsedOptions.lcd = found_lcd; return SSTParsedOptions; } diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index fe12364ff52..f9037b80fe0 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -870,7 +870,15 @@ class CSourcePieceWise_TurbSST final : public CNumerics { ProdDistr[4] = PLim; /*--- Cross diffusion ---*/ - const su2double CrossDiffusionTerm = CDkw_i < 0.0 ? -min(Pw, abs(CrossDiffusionTerm)) : CrossDiffusionTerm; + su2double CrossDiffusionTerm = CDkw_i; + if (sstParsedOptions.lcd) { + su2double exponent = 20.0; + if (sstParsedOptions.version == SST_OPTIONS::V2003) exponent = 10.0; + CrossDiffusionTerm = max(CrossDiffusionTerm, pow(10.0, -exponent)); + } else { + CrossDiffusionTerm = CrossDiffusionTerm < 0.0 ? -min(Pw, abs(CrossDiffusionTerm)) : CrossDiffusionTerm; + } + Residual[1] += (1.0 - F1_i) * CrossDiffusionTerm * Volume; /*--- Contribution due to 2D axisymmetric formulation ---*/ From 9a6fc06707e2bcf0b7d4b15313d40fa9cd4870f7 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 26 Jul 2024 14:48:44 +0200 Subject: [PATCH 19/44] - changed default for cross diffusion --- Common/include/option_structure.hpp | 10 +++++----- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 3ba4fcf3621..2104fe07808 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -994,7 +994,7 @@ enum class SST_OPTIONS { LLT, /*!< \brief Menter k-w SST model with Lower Limits Treatments. */ PRODLIM, /*!< \brief Menter k-w SST model with user-defined production limiter constant. */ NEWBC, /*!< \brief Menter k-w SST model with new boundary conditions. */ - LCD, /*!< \brief Menter k-w SST model with lower limit on the cross-diffusion term. */ + NCD, /*!< \brief Menter k-w SST model with lower limit on the cross-diffusion term. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) @@ -1011,7 +1011,7 @@ static const MapType SST_Options_Map = { MakePair("LLT", SST_OPTIONS::LLT) MakePair("PRODLIM", SST_OPTIONS::PRODLIM) MakePair("NEWBC", SST_OPTIONS::NEWBC) - MakePair("LCD", SST_OPTIONS::LCD) + MakePair("NCD", SST_OPTIONS::NCD) }; /*! @@ -1027,7 +1027,7 @@ struct SST_ParsedOptions { bool llt = false; /*!< \brief Bool for Lower Limits Treatment. */ bool prodLim = false; /*!< \brief Bool for user-defined production limiter constant. */ bool newBC = false; /*!< \brief Bool for new boundary conditions. */ - bool lcd = true; /*!< \brief Bool lower limit on the cross-diffusion term. */ + bool ncd = false; /*!< \brief Bool lower limit on the cross-diffusion term. */ }; /*! @@ -1049,7 +1049,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne const bool found_llt = IsPresent(SST_OPTIONS::LLT); const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM); const bool found_newBC = IsPresent(SST_OPTIONS::NEWBC); - const bool found_lcd = IsPresent(SST_OPTIONS::LCD); + const bool found_ncd = IsPresent(SST_OPTIONS::NCD); const bool found_1994 = IsPresent(SST_OPTIONS::V1994); const bool found_2003 = IsPresent(SST_OPTIONS::V2003); @@ -1096,7 +1096,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.llt = found_llt; SSTParsedOptions.prodLim = found_prodLim; SSTParsedOptions.newBC = found_newBC; - SSTParsedOptions.lcd = found_lcd; + SSTParsedOptions.ncd = found_ncd; return SSTParsedOptions; } diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index f9037b80fe0..8d17dca60c1 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -871,12 +871,12 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Cross diffusion ---*/ su2double CrossDiffusionTerm = CDkw_i; - if (sstParsedOptions.lcd) { + if (sstParsedOptions.ncd) { + CrossDiffusionTerm = CrossDiffusionTerm < 0.0 ? -min(Pw, abs(CrossDiffusionTerm)) : CrossDiffusionTerm; + } else { su2double exponent = 20.0; if (sstParsedOptions.version == SST_OPTIONS::V2003) exponent = 10.0; CrossDiffusionTerm = max(CrossDiffusionTerm, pow(10.0, -exponent)); - } else { - CrossDiffusionTerm = CrossDiffusionTerm < 0.0 ? -min(Pw, abs(CrossDiffusionTerm)) : CrossDiffusionTerm; } Residual[1] += (1.0 - F1_i) * CrossDiffusionTerm * Volume; From 251ae271f51a4c9e4f5780f0838542174bb0bd1c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Aug 2024 18:38:46 +0200 Subject: [PATCH 20/44] - Restored old computation of Cross-Diffusion terms in the source residuals of SST --- Common/include/option_structure.hpp | 5 ----- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 11 +---------- SU2_CFD/src/variables/CTurbSSTVariable.cpp | 6 +++--- 3 files changed, 4 insertions(+), 18 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 2104fe07808..bec6739be30 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -994,7 +994,6 @@ enum class SST_OPTIONS { LLT, /*!< \brief Menter k-w SST model with Lower Limits Treatments. */ PRODLIM, /*!< \brief Menter k-w SST model with user-defined production limiter constant. */ NEWBC, /*!< \brief Menter k-w SST model with new boundary conditions. */ - NCD, /*!< \brief Menter k-w SST model with lower limit on the cross-diffusion term. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) @@ -1011,7 +1010,6 @@ static const MapType SST_Options_Map = { MakePair("LLT", SST_OPTIONS::LLT) MakePair("PRODLIM", SST_OPTIONS::PRODLIM) MakePair("NEWBC", SST_OPTIONS::NEWBC) - MakePair("NCD", SST_OPTIONS::NCD) }; /*! @@ -1027,7 +1025,6 @@ struct SST_ParsedOptions { bool llt = false; /*!< \brief Bool for Lower Limits Treatment. */ bool prodLim = false; /*!< \brief Bool for user-defined production limiter constant. */ bool newBC = false; /*!< \brief Bool for new boundary conditions. */ - bool ncd = false; /*!< \brief Bool lower limit on the cross-diffusion term. */ }; /*! @@ -1049,7 +1046,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne const bool found_llt = IsPresent(SST_OPTIONS::LLT); const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM); const bool found_newBC = IsPresent(SST_OPTIONS::NEWBC); - const bool found_ncd = IsPresent(SST_OPTIONS::NCD); const bool found_1994 = IsPresent(SST_OPTIONS::V1994); const bool found_2003 = IsPresent(SST_OPTIONS::V2003); @@ -1096,7 +1092,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.llt = found_llt; SSTParsedOptions.prodLim = found_prodLim; SSTParsedOptions.newBC = found_newBC; - SSTParsedOptions.ncd = found_ncd; return SSTParsedOptions; } diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 8d17dca60c1..9a028463a06 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -870,16 +870,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { ProdDistr[4] = PLim; /*--- Cross diffusion ---*/ - su2double CrossDiffusionTerm = CDkw_i; - if (sstParsedOptions.ncd) { - CrossDiffusionTerm = CrossDiffusionTerm < 0.0 ? -min(Pw, abs(CrossDiffusionTerm)) : CrossDiffusionTerm; - } else { - su2double exponent = 20.0; - if (sstParsedOptions.version == SST_OPTIONS::V2003) exponent = 10.0; - CrossDiffusionTerm = max(CrossDiffusionTerm, pow(10.0, -exponent)); - } - - Residual[1] += (1.0 - F1_i) * CrossDiffusionTerm * Volume; + Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; /*--- Contribution due to 2D axisymmetric formulation ---*/ diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index de200423ca6..8fa03688303 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -72,14 +72,14 @@ void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_visco // CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -prod_lim_const)); su2double exponent = 20.0; if (sstParsedOptions.version == SST_OPTIONS::V2003) exponent = 10.0; - const su2double CDkw_Clipped = max(CDkw(iPoint), pow(10.0, -exponent)); + CDkw(iPoint)= max(CDkw(iPoint), pow(10.0, -exponent)); /*--- F1 ---*/ arg2A = sqrt(Solution(iPoint,0))/(beta_star*Solution(iPoint,1)*val_dist+EPS*EPS); arg2B = 500.0*val_viscosity / (val_density*val_dist*val_dist*Solution(iPoint,1)+EPS*EPS); arg2 = max(arg2A, arg2B); - arg1 = min(arg2, 4.0*val_density*sigma_om2*Solution(iPoint,0) / (CDkw_Clipped*val_dist*val_dist+EPS*EPS)); + arg1 = min(arg2, 4.0*val_density*sigma_om2*Solution(iPoint,0) / (CDkw(iPoint)*val_dist*val_dist+EPS*EPS)); F1(iPoint) = tanh(pow(arg1, 4.0)); /*--- F2 ---*/ @@ -94,7 +94,7 @@ void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_visco F1(iPoint) = max(F1(iPoint), F3); } - AD::SetPreaccOut(F1(iPoint)); AD::SetPreaccOut(F2(iPoint)); AD::SetPreaccOut(CDkw_Clipped); + AD::SetPreaccOut(F1(iPoint)); AD::SetPreaccOut(F2(iPoint)); AD::SetPreaccOut(CDkw(iPoint)); AD::EndPreacc(); } From 9c7c950a0505c5f592eb038a500ee36c11a100af Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 30 Aug 2024 17:13:40 +0200 Subject: [PATCH 21/44] Added BCs for SST-SUST --- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 3fb642f1dfc..5ba522b22fa 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -133,8 +133,12 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh su2double omega_Inf = rhoInf*kine_Inf/(muLamInf*viscRatio); if (sstParsedOptions.newBC) { - su2double omega_Inf = 10 * sqrt(VelMag2) / config->GetLDomain(); - su2double kine_Inf = omega_Inf*(muLamInf*viscRatio)/rhoInf; + omega_Inf = 10 * sqrt(VelMag2) / config->GetLDomain(); + kine_Inf = omega_Inf*(muLamInf*viscRatio)/rhoInf; + } else if (sstParsedOptions.sust) { + omega_Inf = 5*sqrt(VelMag2)/config->GetLength_Reynolds(); + /*--- not good ---*/ + kine_Inf = pow(10.0, -6) * VelMag2; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ } Solution_Inf[0] = kine_Inf; From 9d82d5198b3c7682baa9d892d9868acd501740f8 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 30 Aug 2024 17:38:34 +0200 Subject: [PATCH 22/44] - Removed duplicate SST options --- Common/include/option_structure.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bf6c0403c26..b8be96692b1 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -995,7 +995,6 @@ enum class SST_OPTIONS { COMP_Sarkar, /*!< \brief Menter k-w SST model with Compressibility correction of Sarkar. */ DLL, /*!< \brief Menter k-w SST model with dimensionless lower limit clipping of turbulence variables. */ FULLPROD, /*!< \brief Menter k-w SST model with full production term. */ - LLT, /*!< \brief Menter k-w SST model with Lower Limits Treatments. */ PRODLIM, /*!< \brief Menter k-w SST model with user-defined production limiter constant. */ NEWBC, /*!< \brief Menter k-w SST model with new boundary conditions. */ }; @@ -1014,7 +1013,6 @@ static const MapType SST_Options_Map = { MakePair("COMPRESSIBILITY-SARKAR", SST_OPTIONS::COMP_Sarkar) MakePair("DIMENSIONLESS_LIMIT", SST_OPTIONS::DLL) MakePair("FULLPROD", SST_OPTIONS::FULLPROD) - MakePair("LLT", SST_OPTIONS::LLT) MakePair("PRODLIM", SST_OPTIONS::PRODLIM) MakePair("NEWBC", SST_OPTIONS::NEWBC) }; @@ -1032,7 +1030,6 @@ struct SST_ParsedOptions { bool compSarkar = false; /*!< \brief Bool for compressibility correction of Sarkar. */ bool dll = false; /*!< \brief Bool dimensionless lower limit. */ bool fullProd = false; /*!< \brief Bool for full production term. */ - bool llt = false; /*!< \brief Bool for Lower Limits Treatment. */ bool prodLim = false; /*!< \brief Bool for user-defined production limiter constant. */ bool newBC = false; /*!< \brief Bool for new boundary conditions. */ }; @@ -1053,7 +1050,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne }; const bool found_fullProd = IsPresent(SST_OPTIONS::FULLPROD); - const bool found_llt = IsPresent(SST_OPTIONS::LLT); const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM); const bool found_newBC = IsPresent(SST_OPTIONS::NEWBC); @@ -1115,7 +1111,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.dll = sst_dll; SSTParsedOptions.fullProd = found_fullProd; - SSTParsedOptions.llt = found_llt; SSTParsedOptions.prodLim = found_prodLim; SSTParsedOptions.newBC = found_newBC; return SSTParsedOptions; From c10b0066d29040b704c85d81fa547bd85be27e74 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 30 Aug 2024 17:51:46 +0200 Subject: [PATCH 23/44] - Removed not used config variables --- Common/include/CConfig.hpp | 4 ---- Common/src/CConfig.cpp | 3 --- 2 files changed, 7 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 06d7882512a..b61b634b169 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1171,8 +1171,6 @@ class CConfig { nHistoryOutput, nVolumeOutput; /*!< \brief Number of variables printed to the history file. */ bool Multizone_Residual; /*!< \brief Determines if memory should be allocated for the multizone residual. */ SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */ - su2double lowerLimitTKE, - lowerLimitDissipation; su2double prodLimConst; su2double LDomain; SA_ParsedOptions saParsedOptions; /*!< \brief Additional parameters for the SA turbulence model. */ @@ -9875,8 +9873,6 @@ class CConfig { */ SST_ParsedOptions GetSSTParsedOptions() const { return sstParsedOptions; } - su2double GetLowerLimitTKE() const { return lowerLimitTKE; } - su2double GetLowerLimitDissipation() const { return lowerLimitDissipation; } su2double GetProdLimConst() const { return prodLimConst; } su2double GetLDomain() const { return LDomain; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index d14a8d758f8..49d9301ece3 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1118,8 +1118,6 @@ void CConfig::SetConfig_Options() { /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SA turbulence model options/corrections. \n Options: see \link SA_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumListOption("SA_OPTIONS", nSA_Options, SA_Options, SA_Options_Map); - addDoubleOption("LOWER_LIMIT_TKE", lowerLimitTKE, 1e-20); - addDoubleOption("LOWER_LIMIT_DISSIPATION", lowerLimitDissipation, 1e-6); addDoubleOption("PROD_LIM_CONST", prodLimConst, 20.0); addDoubleOption("L_DOMAIN", LDomain, 1.0); @@ -6230,7 +6228,6 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "." << endl; if (sstParsedOptions.prodLim) cout << "Changing the value of the TKE production limiter constant to " << prodLimConst << endl; - if (sstParsedOptions.llt) cout << "Changing the value of the lower limits of TKE and Omega " << endl; break; } From be98cdd5f5c14da59c9f126f4b1393cd396d039e Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 30 Aug 2024 18:06:26 +0200 Subject: [PATCH 24/44] - fixed merge error --- SU2_CFD/src/solvers/CEulerSolver.cpp | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index da39dda790d..b1791925cf5 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -7357,6 +7357,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con /*--- Supersonic inlet flow: there are no outgoing characteristics, so all flow variables can be imposed at the inlet. ---*/ + /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -7378,19 +7379,12 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con /*--- Compute the energy from the specified state. ---*/ - const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity); - const su2double Energy_woTKE = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; - su2double Energy = Energy_woTKE; - - /*--- Loop over all the vertices on this boundary marker ---*/ + const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity); + const su2double Energy_woTKE = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; + su2double Energy = Energy_woTKE; - SU2_OMP_FOR_DYN(OMP_MIN_SIZE) - for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { - const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); if (tkeNeeded) Energy = Energy_woTKE + turbNodes->GetSolution(iPoint,0); - if (!geometry->nodes->GetDomain(iPoint)) continue; - /*--- Allocate the value at the inlet ---*/ auto* V_inlet = GetCharacPrimVar(val_marker, iVertex); From fc30f5d2bc1b98d59272e231d2ee72434a7ee377 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 2 Sep 2024 23:34:50 +0200 Subject: [PATCH 25/44] - Include full production and SSTm into the computation of the Stress Tensor --- SU2_CFD/include/numerics/CNumerics.hpp | 10 ++++--- .../include/numerics/flow/flow_diffusion.hpp | 4 ++- .../numerics/turbulent/turb_sources.hpp | 9 ++++--- .../numerics_simd/flow/diffusion/common.hpp | 6 +++-- .../flow/diffusion/viscous_fluxes.hpp | 5 +++- .../include/solvers/CFVMFlowSolverBase.inl | 10 ++++++- SU2_CFD/include/variables/CVariable.hpp | 12 ++++----- .../interfaces/fsi/CFlowTractionInterface.cpp | 10 ++++++- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 26 ++++++++++++------- SU2_CFD/src/numerics/flow/flow_sources.cpp | 5 +++- SU2_CFD/src/solvers/CAdjNSSolver.cpp | 9 ++++++- SU2_CFD/src/solvers/CIncNSSolver.cpp | 4 ++- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 11 +++++++- SU2_CFD/src/solvers/CNSSolver.cpp | 15 ++++++++++- SU2_CFD/src/solvers/CSolver.cpp | 6 ++++- 15 files changed, 107 insertions(+), 35 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 151bcdc2404..050a945b818 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -479,12 +479,14 @@ class CNumerics { template FORCEINLINE static void ComputeStressTensor(size_t nDim, Mat1& stress, const Mat2& velgrad, Scalar viscosity, Scalar density=0.0, - Scalar turb_ke=0.0, bool reynolds3x3=false){ + Scalar turb_ke=0.0, bool fullProd=false, bool modified=false, bool reynolds3x3=false){ Scalar divVel = 0.0; for (size_t iDim = 0; iDim < nDim; iDim++) { divVel += velgrad[iDim][iDim]; } - Scalar pTerm = 2./3. * (divVel * viscosity + density * turb_ke); + Scalar pTerm = 0.0; + if (!modified) pTerm += 2./3. * density * turb_ke; + if (fullProd) pTerm += 2./3. * divVel * viscosity; for (size_t iDim = 0; iDim < nDim; iDim++){ for (size_t jDim = 0; jDim < nDim; jDim++){ @@ -558,9 +560,9 @@ class CNumerics { template NEVERINLINE static void ComputePerturbedRSM(size_t nDim, size_t uq_eigval_comp, bool uq_permute, su2double uq_delta_b, su2double uq_urlx, const Mat1& velgrad, Scalar density, - Scalar viscosity, Scalar turb_ke, Mat2& MeanPerturbedRSM) { + Scalar viscosity, Scalar turb_ke, bool fullProd, bool modified, Mat2& MeanPerturbedRSM) { Scalar MeanReynoldsStress[3][3]; - ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke, true); + ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke, fullProd, modified, true); for (size_t iDim = 0; iDim < 3; iDim++) for (size_t jDim = 0; jDim < 3; jDim++) MeanReynoldsStress[iDim][jDim] /= -density; diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp index a78c23bd270..6612a90fd79 100644 --- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp +++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp @@ -203,7 +203,9 @@ class CAvgGrad_Base : public CNumerics { const su2double* const *val_gradprimvar, su2double val_turb_ke, su2double val_laminar_viscosity, - su2double val_eddy_viscosity); + su2double val_eddy_viscosity, + bool SSTFullProduction, + bool SSTm); /*! * \brief Get a component of the viscous stress tensor. diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 4f55e9f9d3e..662d1d21178 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -837,7 +837,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics { switch (sstParsedOptions.production) { case SST_OPTIONS::UQ: ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, PrimVar_Grad_i + idx.Velocity(), - Density_i, Eddy_Viscosity_i, ScalarVar_i[0], MeanPerturbedRSM); + Density_i, Eddy_Viscosity_i, ScalarVar_i[0], sstParsedOptions.fullProd, + sstParsedOptions.modified, MeanPerturbedRSM); P_Base = PerturbedStrainMag(ScalarVar_i[0]); break; @@ -873,8 +874,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics { const su2double prod_limit = prod_lim_const * beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0]; su2double P = Eddy_Viscosity_i * pow(P_Base, 2); - if (!sstParsedOptions.modified) P -= Eddy_Viscosity_i * diverg*diverg * 2.0/3.0; - if (sstParsedOptions.fullProd) P -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.0; + if (sstParsedOptions.fullProd) P -= Eddy_Viscosity_i * diverg*diverg * 2.0/3.0; + if (!sstParsedOptions.modified) P -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.0; su2double PLim = 0.0; if ( P > prod_limit ) PLim = 1.0; @@ -950,7 +951,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Implicit part ---*/ Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt); - if (sstParsedOptions.fullProd) Jacobian_i[0][0] -= diverg * Volume*2.0/3.0; + if (!sstParsedOptions.modified) Jacobian_i[0][0] -= diverg * Volume*2.0/3.0; Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * (1.0 + zetaFMt); Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume * (1.0 - 0.09/beta_blended * zetaFMt); diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp index 77b4d41fce9..81e1060cb3f 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp @@ -97,7 +97,9 @@ FORCEINLINE MatrixDbl stressTensor(Double viscosity, template NEVERINLINE void addPerturbedRSM(const PrimitiveType& V, const MatrixType& grad, - const Double& turb_ke, + const Double& turb_ke, + const bool FullProd, + const bool modified, MatrixDbl& tau, Ts... uq_args) { /*--- Handle SIMD dimensions 1 by 1. ---*/ @@ -109,7 +111,7 @@ NEVERINLINE void addPerturbedRSM(const PrimitiveType& V, su2double rsm[3][3]; CNumerics::ComputePerturbedRSM(nDim, uq_args..., velgrad, V.density()[k], - V.eddyVisc()[k], turb_ke[k], rsm); + V.eddyVisc()[k], turb_ke[k], FullProd, modified, rsm); for (size_t iDim = 0; iDim < nDim; ++iDim) for (size_t jDim = 0; jDim < nDim; ++jDim) 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 9b34873fcdd..a517d0d2d37 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -140,6 +140,9 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { const auto& solution = static_cast(solution_); const auto& gradient = solution.GetGradient_Primitive(); + const bool SSTm = config.GetSSTParsedOptions().modified; + const bool SSTFullProduction = config.GetSSTParsedOptions().fullProd; + /*--- Compute distance and handle zero without "ifs" by making it large. ---*/ auto dist2_ij = squaredNorm(vector_ij); @@ -158,7 +161,7 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { if(uq) { Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + gatherVariables(jPoint, turbVars->GetSolution())); - addPerturbedRSM(avgV, avgGrad, turb_ke, tau, + addPerturbedRSM(avgV, avgGrad, turb_ke, SSTFullProduction, SSTm, tau, uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx); } diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index f60db438e18..616910bd62c 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2364,6 +2364,9 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr const su2double Prandtl_Lam = config->GetPrandtl_Lam(); const bool energy = config->GetEnergy_Equation(); const bool QCR = config->GetSAParsedOptions().qcr2000; + TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + const bool SSTm = config->GetSSTParsedOptions().modified; + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; const bool axisymmetric = config->GetAxisymmetric(); const bool roughwall = (config->GetnRoughWall() > 0); const bool nemo = config->GetNEMOProblem(); @@ -2447,7 +2450,12 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr } /*--- Evaluate Tau ---*/ - CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity); + su2double tke = 0.0; + // if (TurbFamily == TURB_FAMILY::KW) { + + // } + // Don't knowhow to retrieve the TKE yet + CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity, Density, tke, SSTFullProduction, SSTm); /*--- If necessary evaluate the QCR contribution to Tau ---*/ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 62e7d48a053..d196ab809ad 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1588,6 +1588,12 @@ class CVariable { */ inline virtual su2double GetGradient_Primitive(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return 0.0; } + /*! + * \brief A virtual member. + * \return Value of the primitive variables gradient. + */ + inline virtual CMatrixView GetGradient_Primitive(unsigned long iPoint, unsigned long iVar=0) { return nullptr; } + /*! * \brief Get the primitive variable gradients for all points. * \return Reference to primitive variable gradient. @@ -1609,12 +1615,6 @@ class CVariable { */ inline virtual su2double GetLimiter_Primitive(unsigned long iPoint, unsigned long iVar) const { return 0.0; } - /*! - * \brief A virtual member. - * \return Value of the primitive variables gradient. - */ - inline virtual CMatrixView GetGradient_Primitive(unsigned long iPoint, unsigned long iVar=0) { return nullptr; } - /*! * \brief A virtual member. * \return Value of the primitive variables gradient. diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index f4395640df1..1c542bb0460 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -179,6 +179,9 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry unsigned long Vertex_Flow, unsigned long Point_Struct) { const auto Point_Flow = flow_geometry->vertex[Marker_Flow][Vertex_Flow]->GetNode(); + TURB_FAMILY TurbFamily = TurbModelFamily(flow_config->GetKind_Turb_Model()); + const bool SSTm = flow_config->GetSSTParsedOptions().modified; + const bool SSTFullProduction = flow_config->GetSSTParsedOptions().fullProd; // Get the normal at the vertex: this normal goes inside the fluid domain. const su2double* Normal_Flow = flow_geometry->vertex[Marker_Flow][Vertex_Flow]->GetNormal(); @@ -205,7 +208,12 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow); su2double tau[3][3]; - CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetVelocityGradient(Point_Flow), Viscosity); + su2double TKE = 0.0, Density = 0.0; + if (TurbFamily == TURB_FAMILY::KW && !SSTm) { + TKE = 0.0; // Have to find out how to pass the tke value if needed + Density = 0.0; + } + CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetVelocityGradient(Point_Flow), Viscosity, Density, TKE, SSTFullProduction, SSTm); for (auto iVar = 0u; iVar < nVar; iVar++) { for (auto jVar = 0u; jVar < nVar; jVar++) { // Viscous component in the tn vector --> Units of force (non-dimensional). diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 962445c2851..0780dfaaebb 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -122,7 +122,9 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, const su2double* const *val_gradprimvar, const su2double val_turb_ke, const su2double val_laminar_viscosity, - const su2double val_eddy_viscosity) { + const su2double val_eddy_viscosity, + bool SSTFullProduction, + bool SSTm) { const su2double Density = val_primvar[nDim+2]; @@ -132,7 +134,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, if (sstParsedOptions.uq) { // laminar part - ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity); + ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity, Density, val_turb_ke, SSTFullProduction, SSTm); // add turbulent part which was perturbed for (unsigned short iDim = 0 ; iDim < nDim; iDim++) for (unsigned short jDim = 0 ; jDim < nDim; jDim++) @@ -140,7 +142,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, } else { const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; // turb_ke is not considered in the stress tensor, see #797 - ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, su2double(0.0)); + ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, val_turb_ke, SSTFullProduction, SSTm); } } @@ -370,6 +372,8 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; + const bool SSTm = config->GetSSTParsedOptions().modified; + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; /*--- Normalized normal vector ---*/ @@ -429,13 +433,13 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) if (sstParsedOptions.uq) { ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, Mean_GradPrimVar+1, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, - Mean_turb_ke, MeanPerturbedRSM); + Mean_turb_ke, SSTFullProduction, SSTm, MeanPerturbedRSM); } /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, SSTFullProduction, SSTm); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); @@ -555,6 +559,8 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; + const bool SSTm = config->GetSSTParsedOptions().modified; + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; /*--- Normalized normal vector ---*/ @@ -614,12 +620,12 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi if (sstParsedOptions.uq) { ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, Mean_GradPrimVar+1, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, - Mean_turb_ke, MeanPerturbedRSM); + Mean_turb_ke, SSTFullProduction, SSTm, MeanPerturbedRSM); } /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, SSTFullProduction, SSTm); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); @@ -873,6 +879,8 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; + const bool SSTm = config->GetSSTParsedOptions().modified; + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; /*--- Normalized normal vector ---*/ @@ -944,13 +952,13 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c if (sstParsedOptions.uq) { ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, Mean_GradPrimVar+1, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, - Mean_turb_ke, MeanPerturbedRSM); + Mean_turb_ke, SSTFullProduction, SSTm, MeanPerturbedRSM); } /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, SSTFullProduction, SSTm); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 9ed1305f3ce..9e63716f975 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -253,6 +253,9 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo su2double yinv, Velocity_i[3]; unsigned short iDim, iVar, jVar; + // I don't know how to pass the tke yet + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; + if (Coord_i[1] > EPS) { yinv = 1.0/Coord_i[1]; @@ -317,7 +320,7 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i); /*--- The full stress tensor is needed for variable density ---*/ - ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity); + ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity, 0.0, 0.0, SSTFullProduction); /*--- Viscous terms. ---*/ diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index 3a219b06ea9..bb40e1dbd31 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -592,6 +592,12 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con dp_drv, dp_drw, dp_drE, dH_dr, dH_dru, dH_drv, dH_drw, dH_drE, H, D[3][3], Dd[3], Mach_Inf, eps, scale = 1.0, RefVel2, RefDensity, Mach2Vel, *Velocity_Inf, factor; + TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + const bool SSTm = config->GetSSTParsedOptions().modified; + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; + const auto* turb_solver = solver_container[TURB_SOL]; + const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? turb_solver->GetNodes() : nullptr; + auto *USens = new su2double[nVar]; auto *UnitNormal = new su2double[nDim]; auto *normal_grad_vel = new su2double[nDim]; @@ -762,7 +768,8 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con /*--- Turbulent kinetic energy ---*/ // turb_ke is not considered in the stress tensor, see #797 val_turb_ke = 0.0; - CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke); + if (TurbFamily == TURB_FAMILY::KW && !SSTm) val_turb_ke = Node_Turb->GetSolution(iPoint, 0); + CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke, SSTFullProduction, SSTm); /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/ diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index d8327f1297a..856e0d5eedd 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -632,6 +632,8 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe const su2double kappa = config->GetwallModel_Kappa(); const su2double B = config->GetwallModel_B(); + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (!config->GetViscous_Wall(iMarker)) continue; @@ -703,7 +705,7 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe const su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); su2double Eddy_Visc_Wall = nodes->GetEddyViscosity(iPoint); - CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, 0.0, 0.0, SSTFullProduction); su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim, tau, UnitNormal, TauTangent); diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 5b12013cee2..653f43372c1 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -875,6 +875,13 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, const unsigned short VEL_INDEX = nodes->GetVelIndex(); const unsigned short TVE_INDEX = nodes->GetTveIndex(); + // Is this only at the wall? + TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + const bool SSTm = config->GetSSTParsedOptions().modified; + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; + const auto* turb_solver = solver_container[TURB_SOL]; + const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? turb_solver->GetNodes() : nullptr; + /*--- Loop over boundary points to calculate energy flux ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for(auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -958,7 +965,9 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, /*--- Compute wall shear stress (using the stress tensor) ---*/ su2double Tau[MAXNDIM][MAXNDIM] = {{0.0}}; - CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity); + su2double tke = 0.0; + if (TurbFamily == TURB_FAMILY::KW && !SSTm) tke = Node_Turb->GetSolution(iPoint, 0); + CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity, Density, tke, SSTFullProduction, SSTm); su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim,Tau,UnitNormal,TauTangent); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 29a336c1a26..7d419c635b3 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -332,6 +332,11 @@ void CNSSolver::AddDynamicGridResidualContribution(unsigned long iPoint, unsigne su2double eddy_viscosity = nodes->GetEddyViscosity(iPoint); su2double total_viscosity = laminar_viscosity + eddy_viscosity; + // No config variable here. How to do it? + // TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + // const bool SSTm = config->GetSSTParsedOptions().modified; + // const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; + /*--- Compute the viscous stress tensor ---*/ su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}; @@ -799,6 +804,12 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, const unsigned short max_iter = config->GetwallModel_MaxIter(); const su2double relax = config->GetwallModel_RelFac(); + TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + const bool SSTm = config->GetSSTParsedOptions().modified; + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; + const auto* turb_solver = solver_container[TURB_SOL]; + const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? turb_solver->GetNodes() : nullptr; + /*--- Compute the recovery factor * use Molecular (Laminar) Prandtl number (see Nichols & Nelson, nomenclature ) ---*/ @@ -904,7 +915,9 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, const su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); su2double Eddy_Visc_Wall = nodes->GetEddyViscosity(iPoint); - CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall); + su2double tke = 0.0; + if (TurbFamily == TURB_FAMILY::KW && !SSTm) tke = Node_Turb->GetSolution(iPoint, 0); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, nodes->GetDensity(iPoint), tke, SSTFullProduction, SSTm); su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim, tau, UnitNormal, TauTangent); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index bfbd4edad58..6a417b41bca 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -3989,6 +3989,8 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config) const bool viscous_flow = config->GetViscous(); const su2double Pressure_Inf = config->GetPressure_FreeStreamND(); + TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { @@ -4019,8 +4021,10 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config) // Calculate tn in the fluid nodes for the viscous term if (viscous_flow) { const su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint); + const su2double Density = base_nodes->GetDensity(iPoint); su2double Tau[3][3]; - CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity); + // Here I do not care for modified version of SST since TKE at wall is 0 + CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity, Density, 0.0, SSTFullProduction); for (unsigned short iDim = 0; iDim < nDim; iDim++) { auxForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], Normal); } From a957ec7c775adf4f2680f2b26c3fc8af15705376 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 3 Sep 2024 13:01:16 +0200 Subject: [PATCH 26/44] Updated the freestream values for the print out --- SU2_CFD/src/solvers/CEulerSolver.cpp | 25 ++++++++++++++++----- SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp | 25 ++++++++++++++++----- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 25 ++++++++++++++++----- SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 25 ++++++++++++++++----- 4 files changed, 80 insertions(+), 20 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index b1791925cf5..2a17456f11c 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1039,15 +1039,30 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); Omega_FreeStream = Density_FreeStream*Tke_FreeStream/(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()); + + if (config->GetSSTParsedOptions().newBC) { + Omega_FreeStream = 10 * ModVel_FreeStream / config->GetLDomain(); + Omega_FreeStreamND = 10 * ModVel_FreeStreamND / config->GetLDomain(); // Should it be non-dimensionalized for the Reynolds length? + + Tke_FreeStream = Omega_FreeStream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + Tke_FreeStreamND = Omega_FreeStreamND*(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStreamND; + } else if (config->GetSSTParsedOptions().sust) { + Omega_FreeStream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + Omega_FreeStreamND = 5*ModVel_FreeStreamND / config->GetLength_Reynolds(); // Should it be non-dimensionalized by what? Dimensions are 1/s + + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + Tke_FreeStreamND = pow(10.0, -6) * ModVel_FreeStreamND*ModVel_FreeStreamND; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + config->SetTke_FreeStream(Tke_FreeStream); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + config->SetOmega_FreeStream(Omega_FreeStream); config->SetOmega_FreeStreamND(Omega_FreeStreamND); if (config->GetTurbulenceIntensity_FreeStream() *100 <= 1.3) { diff --git a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp index 73b04c709eb..1788644e3ee 100644 --- a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp +++ b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp @@ -1055,15 +1055,30 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig *config, Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); + + Omega_FreeStream = Density_FreeStream*Tke_FreeStream/(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()); + Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()); + + if (config->GetSSTParsedOptions().newBC) { + Omega_FreeStream = 10 * ModVel_FreeStream / config->GetLDomain(); + Omega_FreeStreamND = 10 * ModVel_FreeStreamND / config->GetLDomain(); // Should it be non-dimensionalized for the Reynolds length? + + Tke_FreeStream = Omega_FreeStream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + Tke_FreeStreamND = Omega_FreeStreamND*(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStreamND; + } else if (config->GetSSTParsedOptions().sust) { + Omega_FreeStream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + Omega_FreeStreamND = 5*ModVel_FreeStreamND / config->GetLength_Reynolds(); // Should it be non-dimensionalized by what? Dimensions are 1/s + + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + Tke_FreeStreamND = pow(10.0, -6) * ModVel_FreeStreamND*ModVel_FreeStreamND; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + config->SetTke_FreeStream(Tke_FreeStream); config->SetTke_FreeStreamND(Tke_FreeStreamND); - Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max((Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()), 1.e-25); config->SetOmega_FreeStream(Omega_FreeStream); - - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max((Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()), 1.e-25); config->SetOmega_FreeStreamND(Omega_FreeStreamND); /*--- Initialize the dimensionless Fluid Model that will be used to solve the dimensionless problem ---*/ diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 5b81394bd9a..1e1e5931c4d 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -448,15 +448,30 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); Omega_FreeStream = Density_FreeStream*Tke_FreeStream/(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()); + + if (config->GetSSTParsedOptions().newBC) { + Omega_FreeStream = 10 * ModVel_FreeStream / config->GetLDomain(); + Omega_FreeStreamND = 10 * ModVel_FreeStreamND / config->GetLDomain(); // Should it be non-dimensionalized for the Reynolds length? + + Tke_FreeStream = Omega_FreeStream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + Tke_FreeStreamND = Omega_FreeStreamND*(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStreamND; + } else if (config->GetSSTParsedOptions().sust) { + Omega_FreeStream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + Omega_FreeStreamND = 5*ModVel_FreeStreamND / config->GetLength_Reynolds(); // Should it be non-dimensionalized by what? Dimensions are 1/s + + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + Tke_FreeStreamND = pow(10.0, -6) * ModVel_FreeStreamND*ModVel_FreeStreamND; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + config->SetTke_FreeStream(Tke_FreeStream); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + config->SetOmega_FreeStream(Omega_FreeStream); config->SetOmega_FreeStreamND(Omega_FreeStreamND); const su2double MassDiffusivityND = config->GetDiffusivity_Constant() / (Velocity_Ref * Length_Ref); diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 776a3baaafd..2ec6f77f2cc 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -1187,15 +1187,30 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); Omega_FreeStream = Density_FreeStream*Tke_FreeStream/(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()); + + if (config->GetSSTParsedOptions().newBC) { + Omega_FreeStream = 10 * ModVel_FreeStream / config->GetLDomain(); + Omega_FreeStreamND = 10 * ModVel_FreeStreamND / config->GetLDomain(); // Should it be non-dimensionalized for the Reynolds length? + + Tke_FreeStream = Omega_FreeStream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + Tke_FreeStreamND = Omega_FreeStreamND*(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStreamND; + } else if (config->GetSSTParsedOptions().sust) { + Omega_FreeStream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + Omega_FreeStreamND = 5*ModVel_FreeStreamND / config->GetLength_Reynolds(); // Should it be non-dimensionalized by what? Dimensions are 1/s + + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + Tke_FreeStreamND = pow(10.0, -6) * ModVel_FreeStreamND*ModVel_FreeStreamND; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + config->SetTke_FreeStream(Tke_FreeStream); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + config->SetOmega_FreeStream(Omega_FreeStream); config->SetOmega_FreeStreamND(Omega_FreeStreamND); /*--- Initialize the dimensionless Fluid Model that will be used to solve the dimensionless problem ---*/ From 554693686d9064f4182fc474dbae73bb511abc8d Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 3 Sep 2024 15:37:52 +0200 Subject: [PATCH 27/44] - Removed unused variables in turb_sources --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 662d1d21178..53ee039c36a 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -881,9 +881,6 @@ class CSourcePieceWise_TurbSST final : public CNumerics { if ( P > prod_limit ) PLim = 1.0; su2double pk = max(0.0, min(P, prod_limit)); - const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; - const su2double zeta = max(ScalarVar_i[1], eddy_visc_var * F2_i / a1); - /*--- Production limiter only for V2003, recompute for V1994. ---*/ su2double pw; if (sstParsedOptions.version == SST_OPTIONS::V1994) { From 1b8ee5bbae72beeadaf3add93258b8d91f2e6c26 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 3 Sep 2024 17:44:56 +0200 Subject: [PATCH 28/44] - Added F2 blending function as output --- SU2_CFD/src/output/CFlowOutput.cpp | 2 ++ SU2_CFD/src/variables/CTurbSSTVariable.cpp | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index a55134a5066..0a50e6b7521 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1259,6 +1259,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ AddVolumeOutput("DESTR_W", "Destr_W", "SOLUTION", "Destruction of rate of dissipation"); AddVolumeOutput("CDkw", "CDkw", "SOLUTION", "Cross-Diffusion term"); AddVolumeOutput("F1", "F1", "SOLUTION", "F1 blending function"); + AddVolumeOutput("F2", "F2", "SOLUTION", "F2 blending function"); break; case TURB_FAMILY::NONE: @@ -1546,6 +1547,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("DESTR_W", iPoint, Node_Turb->GetDestrW(iPoint)); SetVolumeOutputValue("CDkw", iPoint, Node_Turb->GetCrossDiff(iPoint)); SetVolumeOutputValue("F1", iPoint, Node_Turb->GetF1blending(iPoint)); + SetVolumeOutputValue("F2", iPoint, Node_Turb->GetF2blending(iPoint)); SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); if (limiter) { diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index 8fa03688303..378f3ebc246 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -69,7 +69,6 @@ void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_visco for (unsigned long iDim = 0; iDim < nDim; iDim++) CDkw(iPoint) += Gradient(iPoint,0,iDim)*Gradient(iPoint,1,iDim); CDkw(iPoint) *= 2.0*val_density*sigma_om2/Solution(iPoint,1); - // CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -prod_lim_const)); su2double exponent = 20.0; if (sstParsedOptions.version == SST_OPTIONS::V2003) exponent = 10.0; CDkw(iPoint)= max(CDkw(iPoint), pow(10.0, -exponent)); From 6d3586f91956f4217e10fd1a329118972d55f90e Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 3 Sep 2024 17:50:01 +0200 Subject: [PATCH 29/44] - Restore previously removed variables --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 53ee039c36a..662d1d21178 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -881,6 +881,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics { if ( P > prod_limit ) PLim = 1.0; su2double pk = max(0.0, min(P, prod_limit)); + const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; + const su2double zeta = max(ScalarVar_i[1], eddy_visc_var * F2_i / a1); + /*--- Production limiter only for V2003, recompute for V1994. ---*/ su2double pw; if (sstParsedOptions.version == SST_OPTIONS::V1994) { From f27a7cbfdefe7a8c3df1746d3b097bf124c47b7c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 4 Sep 2024 11:28:25 +0200 Subject: [PATCH 30/44] - clean up of variables for Reynolds Stress Tensor computation --- SU2_CFD/include/numerics/CNumerics.hpp | 10 +++--- .../include/numerics/flow/flow_diffusion.hpp | 4 +-- .../numerics/turbulent/turb_sources.hpp | 3 +- .../numerics_simd/flow/diffusion/common.hpp | 4 +-- .../flow/diffusion/viscous_fluxes.hpp | 5 +-- .../include/solvers/CFVMFlowSolverBase.inl | 3 +- .../interfaces/fsi/CFlowTractionInterface.cpp | 3 +- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 32 +++++++++---------- SU2_CFD/src/numerics/flow/flow_sources.cpp | 5 +-- SU2_CFD/src/solvers/CAdjNSSolver.cpp | 6 ++-- SU2_CFD/src/solvers/CIncNSSolver.cpp | 4 +-- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 11 +------ SU2_CFD/src/solvers/CNSSolver.cpp | 6 ++-- 13 files changed, 32 insertions(+), 64 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 050a945b818..151bcdc2404 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -479,14 +479,12 @@ class CNumerics { template FORCEINLINE static void ComputeStressTensor(size_t nDim, Mat1& stress, const Mat2& velgrad, Scalar viscosity, Scalar density=0.0, - Scalar turb_ke=0.0, bool fullProd=false, bool modified=false, bool reynolds3x3=false){ + Scalar turb_ke=0.0, bool reynolds3x3=false){ Scalar divVel = 0.0; for (size_t iDim = 0; iDim < nDim; iDim++) { divVel += velgrad[iDim][iDim]; } - Scalar pTerm = 0.0; - if (!modified) pTerm += 2./3. * density * turb_ke; - if (fullProd) pTerm += 2./3. * divVel * viscosity; + Scalar pTerm = 2./3. * (divVel * viscosity + density * turb_ke); for (size_t iDim = 0; iDim < nDim; iDim++){ for (size_t jDim = 0; jDim < nDim; jDim++){ @@ -560,9 +558,9 @@ class CNumerics { template NEVERINLINE static void ComputePerturbedRSM(size_t nDim, size_t uq_eigval_comp, bool uq_permute, su2double uq_delta_b, su2double uq_urlx, const Mat1& velgrad, Scalar density, - Scalar viscosity, Scalar turb_ke, bool fullProd, bool modified, Mat2& MeanPerturbedRSM) { + Scalar viscosity, Scalar turb_ke, Mat2& MeanPerturbedRSM) { Scalar MeanReynoldsStress[3][3]; - ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke, fullProd, modified, true); + ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke, true); for (size_t iDim = 0; iDim < 3; iDim++) for (size_t jDim = 0; jDim < 3; jDim++) MeanReynoldsStress[iDim][jDim] /= -density; diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp index 6612a90fd79..a78c23bd270 100644 --- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp +++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp @@ -203,9 +203,7 @@ class CAvgGrad_Base : public CNumerics { const su2double* const *val_gradprimvar, su2double val_turb_ke, su2double val_laminar_viscosity, - su2double val_eddy_viscosity, - bool SSTFullProduction, - bool SSTm); + su2double val_eddy_viscosity); /*! * \brief Get a component of the viscous stress tensor. diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 662d1d21178..759af205e3d 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -837,8 +837,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { switch (sstParsedOptions.production) { case SST_OPTIONS::UQ: ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, PrimVar_Grad_i + idx.Velocity(), - Density_i, Eddy_Viscosity_i, ScalarVar_i[0], sstParsedOptions.fullProd, - sstParsedOptions.modified, MeanPerturbedRSM); + Density_i, Eddy_Viscosity_i, ScalarVar_i[0], MeanPerturbedRSM); P_Base = PerturbedStrainMag(ScalarVar_i[0]); break; diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp index 81e1060cb3f..29c3f82a54d 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp @@ -98,8 +98,6 @@ template NEVERINLINE void addPerturbedRSM(const PrimitiveType& V, const MatrixType& grad, const Double& turb_ke, - const bool FullProd, - const bool modified, MatrixDbl& tau, Ts... uq_args) { /*--- Handle SIMD dimensions 1 by 1. ---*/ @@ -111,7 +109,7 @@ NEVERINLINE void addPerturbedRSM(const PrimitiveType& V, su2double rsm[3][3]; CNumerics::ComputePerturbedRSM(nDim, uq_args..., velgrad, V.density()[k], - V.eddyVisc()[k], turb_ke[k], FullProd, modified, rsm); + V.eddyVisc()[k], turb_ke[k], rsm); for (size_t iDim = 0; iDim < nDim; ++iDim) for (size_t jDim = 0; jDim < nDim; ++jDim) 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 a517d0d2d37..9b34873fcdd 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -140,9 +140,6 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { const auto& solution = static_cast(solution_); const auto& gradient = solution.GetGradient_Primitive(); - const bool SSTm = config.GetSSTParsedOptions().modified; - const bool SSTFullProduction = config.GetSSTParsedOptions().fullProd; - /*--- Compute distance and handle zero without "ifs" by making it large. ---*/ auto dist2_ij = squaredNorm(vector_ij); @@ -161,7 +158,7 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { if(uq) { Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + gatherVariables(jPoint, turbVars->GetSolution())); - addPerturbedRSM(avgV, avgGrad, turb_ke, SSTFullProduction, SSTm, tau, + addPerturbedRSM(avgV, avgGrad, turb_ke, tau, uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx); } diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 616910bd62c..d084cd62c0e 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2366,7 +2366,6 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr const bool QCR = config->GetSAParsedOptions().qcr2000; TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); const bool SSTm = config->GetSSTParsedOptions().modified; - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; const bool axisymmetric = config->GetAxisymmetric(); const bool roughwall = (config->GetnRoughWall() > 0); const bool nemo = config->GetNEMOProblem(); @@ -2455,7 +2454,7 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr // } // Don't knowhow to retrieve the TKE yet - CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity, Density, tke, SSTFullProduction, SSTm); + CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity, Density, tke); /*--- If necessary evaluate the QCR contribution to Tau ---*/ diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index 1c542bb0460..45131b653b1 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -181,7 +181,6 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry const auto Point_Flow = flow_geometry->vertex[Marker_Flow][Vertex_Flow]->GetNode(); TURB_FAMILY TurbFamily = TurbModelFamily(flow_config->GetKind_Turb_Model()); const bool SSTm = flow_config->GetSSTParsedOptions().modified; - const bool SSTFullProduction = flow_config->GetSSTParsedOptions().fullProd; // Get the normal at the vertex: this normal goes inside the fluid domain. const su2double* Normal_Flow = flow_geometry->vertex[Marker_Flow][Vertex_Flow]->GetNormal(); @@ -213,7 +212,7 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry TKE = 0.0; // Have to find out how to pass the tke value if needed Density = 0.0; } - CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetVelocityGradient(Point_Flow), Viscosity, Density, TKE, SSTFullProduction, SSTm); + CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetVelocityGradient(Point_Flow), Viscosity, Density, TKE); for (auto iVar = 0u; iVar < nVar; iVar++) { for (auto jVar = 0u; jVar < nVar; jVar++) { // Viscous component in the tn vector --> Units of force (non-dimensional). diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 0780dfaaebb..b1acc6dad9d 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -122,9 +122,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, const su2double* const *val_gradprimvar, const su2double val_turb_ke, const su2double val_laminar_viscosity, - const su2double val_eddy_viscosity, - bool SSTFullProduction, - bool SSTm) { + const su2double val_eddy_viscosity) { const su2double Density = val_primvar[nDim+2]; @@ -134,7 +132,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, if (sstParsedOptions.uq) { // laminar part - ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity, Density, val_turb_ke, SSTFullProduction, SSTm); + ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity, Density, val_turb_ke); // add turbulent part which was perturbed for (unsigned short iDim = 0 ; iDim < nDim; iDim++) for (unsigned short jDim = 0 ; jDim < nDim; jDim++) @@ -142,7 +140,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, } else { const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; // turb_ke is not considered in the stress tensor, see #797 - ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, val_turb_ke, SSTFullProduction, SSTm); + ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, val_turb_ke); } } @@ -373,7 +371,6 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) unsigned short iVar, jVar, iDim; const bool SSTm = config->GetSSTParsedOptions().modified; - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; /*--- Normalized normal vector ---*/ @@ -407,6 +404,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) Mean_Laminar_Viscosity = 0.5*(Laminar_Viscosity_i + Laminar_Viscosity_j); Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); + const su2double Mean_turb_ke_ForStressTensor = !SSTm ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ @@ -433,13 +431,13 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) if (sstParsedOptions.uq) { ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, Mean_GradPrimVar+1, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, - Mean_turb_ke, SSTFullProduction, SSTm, MeanPerturbedRSM); + Mean_turb_ke, MeanPerturbedRSM); } /*--- Get projected flux tensor (viscous residual) ---*/ - SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, SSTFullProduction, SSTm); + SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke_ForStressTensor, + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); @@ -560,7 +558,6 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi unsigned short iVar, jVar, iDim; const bool SSTm = config->GetSSTParsedOptions().modified; - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; /*--- Normalized normal vector ---*/ @@ -596,6 +593,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); Mean_Thermal_Conductivity = 0.5*(Thermal_Conductivity_i + Thermal_Conductivity_j); + const su2double Mean_turb_ke_ForStressTensor = !SSTm ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ @@ -620,12 +618,12 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi if (sstParsedOptions.uq) { ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, Mean_GradPrimVar+1, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, - Mean_turb_ke, SSTFullProduction, SSTm, MeanPerturbedRSM); + Mean_turb_ke, MeanPerturbedRSM); } /*--- Get projected flux tensor (viscous residual) ---*/ - SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, SSTFullProduction, SSTm); + SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke_ForStressTensor, + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); @@ -880,7 +878,6 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c unsigned short iVar, jVar, iDim; const bool SSTm = config->GetSSTParsedOptions().modified; - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; /*--- Normalized normal vector ---*/ @@ -926,6 +923,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); Mean_Thermal_Conductivity = 0.5*(Thermal_Conductivity_i + Thermal_Conductivity_j); Mean_Cp = 0.5*(Cp_i + Cp_j); + const su2double Mean_turb_ke_ForStressTensor = !SSTm ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ @@ -952,13 +950,13 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c if (sstParsedOptions.uq) { ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, Mean_GradPrimVar+1, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, - Mean_turb_ke, SSTFullProduction, SSTm, MeanPerturbedRSM); + Mean_turb_ke, MeanPerturbedRSM); } /*--- Get projected flux tensor (viscous residual) ---*/ - SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, SSTFullProduction, SSTm); + SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke_ForStressTensor, + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 9e63716f975..51908bd814f 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -253,9 +253,6 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo su2double yinv, Velocity_i[3]; unsigned short iDim, iVar, jVar; - // I don't know how to pass the tke yet - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; - if (Coord_i[1] > EPS) { yinv = 1.0/Coord_i[1]; @@ -320,7 +317,7 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i); /*--- The full stress tensor is needed for variable density ---*/ - ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity, 0.0, 0.0, SSTFullProduction); + ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity, 0.0, 0.0); /*--- Viscous terms. ---*/ diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index bb40e1dbd31..8068cd50bd1 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -594,9 +594,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); const bool SSTm = config->GetSSTParsedOptions().modified; - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; - const auto* turb_solver = solver_container[TURB_SOL]; - const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? turb_solver->GetNodes() : nullptr; + const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? solver_container[TURB_SOL]->GetNodes() : nullptr; auto *USens = new su2double[nVar]; auto *UnitNormal = new su2double[nDim]; @@ -769,7 +767,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con // turb_ke is not considered in the stress tensor, see #797 val_turb_ke = 0.0; if (TurbFamily == TURB_FAMILY::KW && !SSTm) val_turb_ke = Node_Turb->GetSolution(iPoint, 0); - CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke, SSTFullProduction, SSTm); + CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke); /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/ diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 856e0d5eedd..006c8fff098 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -632,8 +632,6 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe const su2double kappa = config->GetwallModel_Kappa(); const su2double B = config->GetwallModel_B(); - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; - for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (!config->GetViscous_Wall(iMarker)) continue; @@ -705,7 +703,7 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe const su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); su2double Eddy_Visc_Wall = nodes->GetEddyViscosity(iPoint); - CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, 0.0, 0.0, SSTFullProduction); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, 0.0); su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim, tau, UnitNormal, TauTangent); diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 653f43372c1..d22e0f956ac 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -875,13 +875,6 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, const unsigned short VEL_INDEX = nodes->GetVelIndex(); const unsigned short TVE_INDEX = nodes->GetTveIndex(); - // Is this only at the wall? - TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); - const bool SSTm = config->GetSSTParsedOptions().modified; - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; - const auto* turb_solver = solver_container[TURB_SOL]; - const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? turb_solver->GetNodes() : nullptr; - /*--- Loop over boundary points to calculate energy flux ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for(auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -965,9 +958,7 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, /*--- Compute wall shear stress (using the stress tensor) ---*/ su2double Tau[MAXNDIM][MAXNDIM] = {{0.0}}; - su2double tke = 0.0; - if (TurbFamily == TURB_FAMILY::KW && !SSTm) tke = Node_Turb->GetSolution(iPoint, 0); - CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity, Density, tke, SSTFullProduction, SSTm); + CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity, Density); // no need for TKE as it is on the wall su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim,Tau,UnitNormal,TauTangent); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 7d419c635b3..b8eec754f4d 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -806,9 +806,7 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); const bool SSTm = config->GetSSTParsedOptions().modified; - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; - const auto* turb_solver = solver_container[TURB_SOL]; - const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? turb_solver->GetNodes() : nullptr; + const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? solver_container[TURB_SOL]->GetNodes() : nullptr; /*--- Compute the recovery factor * use Molecular (Laminar) Prandtl number (see Nichols & Nelson, nomenclature ) ---*/ @@ -917,7 +915,7 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, su2double tke = 0.0; if (TurbFamily == TURB_FAMILY::KW && !SSTm) tke = Node_Turb->GetSolution(iPoint, 0); - CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, nodes->GetDensity(iPoint), tke, SSTFullProduction, SSTm); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, nodes->GetDensity(iPoint), tke); su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim, tau, UnitNormal, TauTangent); From 448a532387b1291815cea6320f43af6e54ee3af5 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 4 Sep 2024 19:28:39 +0200 Subject: [PATCH 31/44] - Use full tke production term in Pw instead of only P_base --- 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 759af205e3d..be2615c5e2f 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -886,7 +886,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Production limiter only for V2003, recompute for V1994. ---*/ su2double pw; if (sstParsedOptions.version == SST_OPTIONS::V1994) { - pw = alfa_blended * Density_i * pow(P_Base, 2); + pw = alfa_blended * Density_i * P / Eddy_Viscosity_i; } else { pw = (alfa_blended * Density_i / Eddy_Viscosity_i) * pk; } From e9fa4e7c28b7f5dba0563d590e91284ff633b932 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 6 Sep 2024 18:11:23 +0200 Subject: [PATCH 32/44] - Fixed Errors in AD compiling --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 2 +- SU2_CFD/src/solvers/CIncNSSolver.cpp | 2 +- SU2_CFD/src/solvers/CSolver.cpp | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 51908bd814f..9ed1305f3ce 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -317,7 +317,7 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i); /*--- The full stress tensor is needed for variable density ---*/ - ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity, 0.0, 0.0); + ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity); /*--- Viscous terms. ---*/ diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 006c8fff098..d8327f1297a 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -703,7 +703,7 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe const su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); su2double Eddy_Visc_Wall = nodes->GetEddyViscosity(iPoint); - CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, 0.0); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall); su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim, tau, UnitNormal, TauTangent); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 6a417b41bca..6bd045c62b3 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4024,7 +4024,8 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config) const su2double Density = base_nodes->GetDensity(iPoint); su2double Tau[3][3]; // Here I do not care for modified version of SST since TKE at wall is 0 - CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity, Density, 0.0, SSTFullProduction); + su2double tke = 0.0; + CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity, Density, tke, SSTFullProduction); for (unsigned short iDim = 0; iDim < nDim; iDim++) { auxForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], Normal); } From c2462e80bb661bec442abd65ce2758df5244265b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 9 Sep 2024 11:48:02 +0200 Subject: [PATCH 33/44] - added tke to the numerics simd computations - TODO: change the jacobian of the stress tensor --- SU2_CFD/include/numerics/CNumerics.hpp | 2 +- .../include/numerics_simd/flow/diffusion/common.hpp | 5 +++-- .../numerics_simd/flow/diffusion/viscous_fluxes.hpp | 11 +++++++---- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 4 ++-- SU2_CFD/src/solvers/CNSSolver.cpp | 9 +++++---- SU2_CFD/src/solvers/CSolver.cpp | 4 +--- 6 files changed, 19 insertions(+), 16 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 151bcdc2404..55d214698c6 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -560,7 +560,7 @@ class CNumerics { su2double uq_urlx, const Mat1& velgrad, Scalar density, Scalar viscosity, Scalar turb_ke, Mat2& MeanPerturbedRSM) { Scalar MeanReynoldsStress[3][3]; - ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke, true); + ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke); for (size_t iDim = 0; iDim < 3; iDim++) for (size_t jDim = 0; jDim < 3; jDim++) MeanReynoldsStress[iDim][jDim] /= -density; diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp index 29c3f82a54d..d4da1f61cf6 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp @@ -71,13 +71,13 @@ FORCEINLINE void correctGradient(const PrimitiveType& V, */ template FORCEINLINE MatrixDbl stressTensor(Double viscosity, - const MatrixDbl& grad) { + const MatrixDbl& grad, Double density=0.0, Double tke=0.0) { /*--- Hydrostatic term. ---*/ Double velDiv = 0.0; for (size_t iDim = 0; iDim < nDim; ++iDim) { velDiv += grad(iDim+1,iDim); } - Double pTerm = 2.0/3.0 * viscosity * velDiv; + Double pTerm = 2.0/3.0 * (viscosity * velDiv + density * tke); MatrixDbl tau; for (size_t iDim = 0; iDim < nDim; ++iDim) { @@ -186,6 +186,7 @@ template FORCEINLINE MatrixDbl stressTensorJacobian(const PrimitiveType& V, const VectorDbl& normal, Double dist_ij) { + /// TODO: include tke if the SST model is used. Double viscosity = V.laminarVisc() + V.eddyVisc(); Double xi = viscosity / (V.density() * dist_ij); MatrixDbl jac; 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 9b34873fcdd..de32485d1ed 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -140,6 +140,8 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { const auto& solution = static_cast(solution_); const auto& gradient = solution.GetGradient_Primitive(); + const bool tkeNeeded = (config.GetKind_Turb_Model() == TURB_MODEL::SST) && !(config.GetSSTParsedOptions().modified); + /*--- Compute distance and handle zero without "ifs" by making it large. ---*/ auto dist2_ij = squaredNorm(vector_ij); @@ -152,12 +154,13 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { if(correct) correctGradient(V, vector_ij, dist2_ij, avgGrad); /*--- Stress and heat flux tensors. ---*/ - - auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad); + Double turb_ke = 0.0; + Double density = avgV.density(); + if (tkeNeeded) turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + + gatherVariables(jPoint, turbVars->GetSolution())); + auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad, density, turb_ke); if(useSA_QCR) addQCR(avgGrad, tau); if(uq) { - Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + - gatherVariables(jPoint, turbVars->GetSolution())); addPerturbedRSM(avgV, avgGrad, turb_ke, tau, uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx); } diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index b1acc6dad9d..8c3984b2f51 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -370,7 +370,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; - const bool SSTm = config->GetSSTParsedOptions().modified; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified); /*--- Normalized normal vector ---*/ @@ -404,7 +404,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) Mean_Laminar_Viscosity = 0.5*(Laminar_Viscosity_i + Laminar_Viscosity_j); Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); - const su2double Mean_turb_ke_ForStressTensor = !SSTm ? Mean_turb_ke : 0.0; + const su2double Mean_turb_ke_ForStressTensor = tkeNeeded ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index b8eec754f4d..fda3e1af390 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -340,7 +340,8 @@ void CNSSolver::AddDynamicGridResidualContribution(unsigned long iPoint, unsigne /*--- Compute the viscous stress tensor ---*/ su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}; - CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), total_viscosity); + su2double tke=0.0, density=0.0; + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), total_viscosity, density, tke); /*--- Dot product of the stress tensor with the grid velocity ---*/ @@ -804,9 +805,9 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, const unsigned short max_iter = config->GetwallModel_MaxIter(); const su2double relax = config->GetwallModel_RelFac(); - TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); const bool SSTm = config->GetSSTParsedOptions().modified; - const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? solver_container[TURB_SOL]->GetNodes() : nullptr; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !SSTm; + const auto* Node_Turb = (tkeNeeded) ? solver_container[TURB_SOL]->GetNodes() : nullptr; /*--- Compute the recovery factor * use Molecular (Laminar) Prandtl number (see Nichols & Nelson, nomenclature ) ---*/ @@ -914,7 +915,7 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, su2double Eddy_Visc_Wall = nodes->GetEddyViscosity(iPoint); su2double tke = 0.0; - if (TurbFamily == TURB_FAMILY::KW && !SSTm) tke = Node_Turb->GetSolution(iPoint, 0); + if (tkeNeeded) tke = Node_Turb->GetSolution(iPoint, 0); CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, nodes->GetDensity(iPoint), tke); su2double TauTangent[MAXNDIM] = {0.0}; diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 6bd045c62b3..25da1e0155a 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4021,11 +4021,9 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config) // Calculate tn in the fluid nodes for the viscous term if (viscous_flow) { const su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint); - const su2double Density = base_nodes->GetDensity(iPoint); su2double Tau[3][3]; // Here I do not care for modified version of SST since TKE at wall is 0 - su2double tke = 0.0; - CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity, Density, tke, SSTFullProduction); + CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity); for (unsigned short iDim = 0; iDim < nDim; iDim++) { auxForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], Normal); } From 51112feb9604217fa9cb1f329d5552c2a7578a3c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 9 Sep 2024 13:54:09 +0200 Subject: [PATCH 34/44] - Fixed UQ problem --- SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp | 2 ++ 1 file changed, 2 insertions(+) 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 de32485d1ed..667f1a6fa61 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -161,6 +161,8 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad, density, turb_ke); if(useSA_QCR) addQCR(avgGrad, tau); if(uq) { + turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + + gatherVariables(jPoint, turbVars->GetSolution())); addPerturbedRSM(avgV, avgGrad, turb_ke, tau, uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx); } From 576e8ee9c4e0b28acaf7b731df099b558fe93bde Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 9 Sep 2024 18:41:00 +0200 Subject: [PATCH 35/44] - fix UQ implementation --- SU2_CFD/include/numerics/CNumerics.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 55d214698c6..151bcdc2404 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -560,7 +560,7 @@ class CNumerics { su2double uq_urlx, const Mat1& velgrad, Scalar density, Scalar viscosity, Scalar turb_ke, Mat2& MeanPerturbedRSM) { Scalar MeanReynoldsStress[3][3]; - ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke); + ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke, true); for (size_t iDim = 0; iDim < 3; iDim++) for (size_t jDim = 0; jDim < 3; jDim++) MeanReynoldsStress[iDim][jDim] /= -density; From acc21ec0ae4f40d6b4648b065a62250997442894 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 9 Sep 2024 19:27:34 +0200 Subject: [PATCH 36/44] - Removed unused variables --- SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp | 2 +- SU2_CFD/src/solvers/CAdjNSSolver.cpp | 6 +++--- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 2 +- SU2_CFD/src/solvers/CNSSolver.cpp | 8 +------- SU2_CFD/src/solvers/CSolver.cpp | 2 -- 5 files changed, 6 insertions(+), 14 deletions(-) diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp index d4da1f61cf6..cabb7f5c793 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp @@ -97,7 +97,7 @@ FORCEINLINE MatrixDbl stressTensor(Double viscosity, template NEVERINLINE void addPerturbedRSM(const PrimitiveType& V, const MatrixType& grad, - const Double& turb_ke, + const Double& turb_ke, MatrixDbl& tau, Ts... uq_args) { /*--- Handle SIMD dimensions 1 by 1. ---*/ diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index 8068cd50bd1..8cf8c23f1da 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -592,9 +592,9 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con dp_drv, dp_drw, dp_drE, dH_dr, dH_dru, dH_drv, dH_drw, dH_drE, H, D[3][3], Dd[3], Mach_Inf, eps, scale = 1.0, RefVel2, RefDensity, Mach2Vel, *Velocity_Inf, factor; - TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); const bool SSTm = config->GetSSTParsedOptions().modified; - const auto* Node_Turb = (TurbFamily == TURB_FAMILY::KW) ? solver_container[TURB_SOL]->GetNodes() : nullptr; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !SSTm; + const auto* Node_Turb = (tkeNeeded) ? solver_container[TURB_SOL]->GetNodes() : nullptr; auto *USens = new su2double[nVar]; auto *UnitNormal = new su2double[nDim]; @@ -766,7 +766,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con /*--- Turbulent kinetic energy ---*/ // turb_ke is not considered in the stress tensor, see #797 val_turb_ke = 0.0; - if (TurbFamily == TURB_FAMILY::KW && !SSTm) val_turb_ke = Node_Turb->GetSolution(iPoint, 0); + if (tkeNeeded) val_turb_ke = Node_Turb->GetSolution(iPoint, 0); CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke); /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/ diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index d22e0f956ac..2183ae32fb1 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -958,7 +958,7 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, /*--- Compute wall shear stress (using the stress tensor) ---*/ su2double Tau[MAXNDIM][MAXNDIM] = {{0.0}}; - CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity, Density); // no need for TKE as it is on the wall + CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity); // no need for TKE as it is on the wall su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim,Tau,UnitNormal,TauTangent); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index fda3e1af390..f2f1d6f4c4e 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -332,16 +332,10 @@ void CNSSolver::AddDynamicGridResidualContribution(unsigned long iPoint, unsigne su2double eddy_viscosity = nodes->GetEddyViscosity(iPoint); su2double total_viscosity = laminar_viscosity + eddy_viscosity; - // No config variable here. How to do it? - // TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); - // const bool SSTm = config->GetSSTParsedOptions().modified; - // const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; - /*--- Compute the viscous stress tensor ---*/ su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}; - su2double tke=0.0, density=0.0; - CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), total_viscosity, density, tke); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), total_viscosity); /*--- Dot product of the stress tensor with the grid velocity ---*/ diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 25da1e0155a..227692d0b83 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -3989,8 +3989,6 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config) const bool viscous_flow = config->GetViscous(); const su2double Pressure_Inf = config->GetPressure_FreeStreamND(); - TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); - const bool SSTFullProduction = config->GetSSTParsedOptions().fullProd; for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { From c6376cd2f671fd88c9db12f1bf94ce2603ed6ed9 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 12 Sep 2024 17:31:21 +0200 Subject: [PATCH 37/44] - Started including tke when computing the speed of sound --- .../flow/convection/centered.hpp | 11 ++++++++++ .../numerics_simd/flow/convection/common.hpp | 15 +++++++++++-- .../numerics_simd/flow/convection/roe.hpp | 21 +++++++++++++++++-- .../numerics_simd/flow/diffusion/common.hpp | 1 - .../flow/diffusion/viscous_fluxes.hpp | 4 ++++ .../include/numerics_simd/flow/variables.hpp | 14 ++++++++++++- 6 files changed, 60 insertions(+), 6 deletions(-) diff --git a/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp b/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp index ae5c23c259c..9b37f976d9a 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp @@ -53,6 +53,8 @@ class CCenteredBase : public Base { const bool dynamicGrid; const su2double stretchParam = 0.3; + using Base::turbVars; + /*! * \brief Constructor, store some constants and forward args to base. */ @@ -96,6 +98,8 @@ class CCenteredBase : public Base { const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT); const auto& solution = static_cast(solution_); + const bool tkeNeeded = config.GetKind_Turb_Model() == TURB_MODEL::SST; + const auto iPoint = geometry.edges->GetNode(iEdge,0); const auto jPoint = geometry.edges->GetNode(iEdge,1); @@ -119,6 +123,13 @@ class CCenteredBase : public Base { avgV.all(iVar) = 0.5 * (V.i.all(iVar) + V.j.all(iVar)); } + if (tkeNeeded) { + V.i.allTurb = gatherVariables<1>(iPoint, turbVars->GetSolution()); + V.j.allTurb = gatherVariables<1>(jPoint, turbVars->GetSolution()); + + avgV.allTurb(0) = 0.5*(V.i.allTurb(0)+V.j.allTurb(0)); + } + /*--- Compute conservative variables. ---*/ CPair > U; diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index e1807f99f6b..e3781a0c09d 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -103,20 +103,30 @@ FORCEINLINE void musclEdgeLimited(Int iPoint, template FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int jPoint, bool muscl, LIMITER limiterType, + bool musclTurb, LIMITER limiterTypeTurb, const CPair& V1st, const VectorDbl& vector_ij, - const VariableType& solution) { + const VariableType& solution, + const bool tkeNeeded) { static_assert(ReconVarType::nVar <= PrimVarType::nVar,""); const auto& gradients = solution.GetGradient_Reconstruction(); const auto& limiters = solution.GetLimiter_Primitive(); + // const auto& gradientsTurb = turbSolution.GetGradient_Reconstruction(); + // const auto& limitersTurb = turbSolution.GetLimiter_Primitive(); + CPair V; for (size_t iVar = 0; iVar < ReconVarType::nVar; ++iVar) { V.i.all(iVar) = V1st.i.all(iVar); V.j.all(iVar) = V1st.j.all(iVar); } + // Only first order for turbulence + if (tkeNeeded) { + V.i.allTurb(0) = V1st.i.allTurb(0); + V.j.allTurb(0) = V1st.j.allTurb(0); + } if (muscl) { switch (limiterType) { @@ -143,9 +153,10 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int for (size_t iDim = 0; iDim < nDim; ++iDim) { v_squared += pow(R*V.j.velocity(iDim) + V.i.velocity(iDim), 2); } + Double tke = R*V1st.j.allTurb(0) + V1st.i.allTurb(0); /*--- Multiply enthalpy by R+1 since v^2 was not divided by (R+1)^2. * Note: a = sqrt((gamma-1) * (H - 0.5 * v^2)) ---*/ - const Double neg_sound_speed = enthalpy * (R+1) < 0.5 * v_squared; + const Double neg_sound_speed = enthalpy * (R+1) < (0.5 * v_squared - tke); /*--- Revert to first order if the state is non-physical. ---*/ Double bad_recon = fmax(neg_p_or_rho, neg_sound_speed); diff --git a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp index d128ae4cbca..dabb24837cf 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp @@ -60,7 +60,11 @@ class CRoeBase : public Base { const bool finestGrid; const bool dynamicGrid; const bool muscl; + const bool musclTurb; const LIMITER typeLimiter; + const LIMITER typeLimiterTurb; + + using Base::turbVars; /*! * \brief Constructor, store some constants and forward args to base. @@ -73,7 +77,9 @@ class CRoeBase : public Base { finestGrid(iMesh == MESH_0), dynamicGrid(config.GetDynamic_Grid()), muscl(finestGrid && config.GetMUSCL_Flow()), - typeLimiter(config.GetKind_SlopeLimit_Flow()) { + musclTurb(finestGrid && config.GetMUSCL_Turb()), + typeLimiter(config.GetKind_SlopeLimit_Flow()), + typeLimiterTurb(config.GetKind_SlopeLimit_Turb()) { } public: @@ -96,9 +102,13 @@ class CRoeBase : public Base { const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT); const auto& solution = static_cast(solution_); + const bool tkeNeeded = config.GetKind_Turb_Model() == TURB_MODEL::SST; + const auto iPoint = geometry.edges->GetNode(iEdge,0); const auto jPoint = geometry.edges->GetNode(iEdge,1); + // cout << gatherVariables(iPoint, turbVars->GetSolution()) << endl; + /*--- Geometric properties. ---*/ const auto vector_ij = distanceVector(iPoint, jPoint, geometry.nodes->GetCoord()); @@ -116,8 +126,13 @@ class CRoeBase : public Base { V1st.i.all = gatherVariables(iPoint, solution.GetPrimitive()); V1st.j.all = gatherVariables(jPoint, solution.GetPrimitive()); + if (tkeNeeded) { + V1st.i.allTurb = gatherVariables(iPoint, turbVars->GetSolution()); + V1st.j.allTurb = gatherVariables(jPoint, turbVars->GetSolution()); + } + auto V = reconstructPrimitives >( - iEdge, iPoint, jPoint, muscl, typeLimiter, V1st, vector_ij, solution); + iEdge, iPoint, jPoint, muscl, typeLimiter, musclTurb, typeLimiterTurb, V1st, vector_ij, solution, tkeNeeded); /*--- Compute conservative variables. ---*/ @@ -229,6 +244,8 @@ class CRoeScheme : public CRoeBase,Decorator> { using Base::kappa; const ENUM_ROELOWDISS typeDissip; + using Base::turbVars; + public: /*! * \brief Constructor, store some constants and forward to base. diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp index cabb7f5c793..5555ca89d2e 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp @@ -186,7 +186,6 @@ template FORCEINLINE MatrixDbl stressTensorJacobian(const PrimitiveType& V, const VectorDbl& normal, Double dist_ij) { - /// TODO: include tke if the SST model is used. Double viscosity = V.laminarVisc() + V.eddyVisc(); Double xi = viscosity / (V.density() * dist_ij); MatrixDbl jac; 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 667f1a6fa61..d20edd4eefc 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -50,6 +50,7 @@ class CNoViscousFlux : public CNumericsSIMD { protected: static constexpr size_t nDim = NDIM; static constexpr size_t nPrimVar = 0; + const CVariable* turbVars = nullptr; template CNoViscousFlux(Ts&...) {} @@ -268,6 +269,7 @@ class CCompressibleViscousFlux : public CCompressibleViscousFluxBase >; using Base::prandtlTurb; + using Base::turbVars; +; /*! * \brief Constructor, initialize constants and booleans. diff --git a/SU2_CFD/include/numerics_simd/flow/variables.hpp b/SU2_CFD/include/numerics_simd/flow/variables.hpp index 17e6f4a090f..44fd815911d 100644 --- a/SU2_CFD/include/numerics_simd/flow/variables.hpp +++ b/SU2_CFD/include/numerics_simd/flow/variables.hpp @@ -50,12 +50,16 @@ struct CCompressiblePrimitives { FORCEINLINE const Double& velocity(size_t iDim) const { return all(iDim+1); } FORCEINLINE const Double* velocity() const { return &velocity(0); } + VectorDbl<1> allTurb; + FORCEINLINE const Double& tke() const { return allTurb(0); } + /*--- Un-reconstructed variables. ---*/ FORCEINLINE Double& speedSound() { return all(nDim+4); } FORCEINLINE Double& laminarVisc() { return all(nDim+5); } FORCEINLINE Double& eddyVisc() { return all(nDim+6); } FORCEINLINE Double& thermalCond() { return all(nDim+7); } FORCEINLINE Double& cp() { return all(nDim+8); } + FORCEINLINE Double& tke() { return all(nDim+9); } FORCEINLINE const Double& speedSound() const { return all(nDim+4); } FORCEINLINE const Double& laminarVisc() const { return all(nDim+5); } FORCEINLINE const Double& eddyVisc() const { return all(nDim+6); } @@ -70,17 +74,22 @@ template struct CCompressibleConservatives { static constexpr size_t nDim = nDim_; static constexpr size_t nVar = nDim+2; + static constexpr size_t nVarTurb = 1; VectorDbl all; + VectorDbl allTurb; FORCEINLINE Double& density() { return all(0); } FORCEINLINE Double& rhoEnergy() { return all(nDim+1); } FORCEINLINE Double& momentum(size_t iDim) { return all(iDim+1); } + FORCEINLINE Double& rhoTke() { return allTurb(0); } FORCEINLINE const Double& density() const { return all(0); } FORCEINLINE const Double& rhoEnergy() const { return all(nDim+1); } FORCEINLINE const Double& momentum(size_t iDim) const { return all(iDim+1); } + FORCEINLINE const Double& rhoTke() const { return allTurb(0); } FORCEINLINE Double energy() const { return rhoEnergy() / density(); } FORCEINLINE const Double* momentum() const { return &momentum(0); } + FORCEINLINE Double tke() const { return rhoTke() / density(); } }; /*! @@ -94,6 +103,7 @@ FORCEINLINE CCompressibleConservatives compressibleConservatives(const CCo U.momentum(iDim) = V.density() * V.velocity(iDim); } U.rhoEnergy() = V.density() * V.enthalpy() - V.pressure(); + U.rhoTke() = V.density() * V.tke(); return U; } @@ -105,6 +115,7 @@ struct CRoeVariables { Double density; VectorDbl velocity; Double enthalpy; + Double tke; Double speedSound; Double projVel; }; @@ -124,7 +135,8 @@ FORCEINLINE CRoeVariables roeAveragedVariables(Double gamma, roeAvg.velocity(iDim) = (R*V.j.velocity(iDim) + V.i.velocity(iDim)) * D; } roeAvg.enthalpy = (R*V.j.enthalpy() + V.i.enthalpy()) * D; - roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity))); + roeAvg.tke = (R*V.j.tke() + V.i.tke()) * D; + roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity) - roeAvg.tke)); roeAvg.projVel = dot(roeAvg.velocity, normal); return roeAvg; } From c5aeb5af13797b31cf87df17307787c4a0e24b60 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 16 Sep 2024 20:42:15 +0200 Subject: [PATCH 38/44] - Start to include tke only when not m version of SST is used --- Common/include/CConfig.hpp | 2 +- .../include/solvers/CFVMFlowSolverBase.inl | 2 +- SU2_CFD/src/solvers/CEulerSolver.cpp | 50 +++++++++++++------ SU2_CFD/src/solvers/CNSSolver.cpp | 5 +- 4 files changed, 39 insertions(+), 20 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 7a46555c8a6..deac682f12f 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1961,7 +1961,7 @@ class CConfig { * \brief Get the value of the non-dimensionalized freestream energy. * \return Non-dimensionalized freestream energy. */ - su2double GetEnergy_FreeStreamND(void) const { return Energy_FreeStreamND; } + su2double GetEnergy_FreeStreamND(void) const { cout << "La chiedo non-dimensionale " << endl; return Energy_FreeStreamND; } /*! * \brief Get the value of the non-dimensionalized freestream viscosity. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index d084cd62c0e..da5f2cc7e78 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -432,7 +432,7 @@ void CFVMFlowSolverBase::Viscous_Residual_impl(unsigned long iEdge, CGeome CNumerics *numerics, CConfig *config) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); CVariable* turbNodes = nullptr; if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 0069db73130..1a2b494eea1 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -64,6 +64,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING); const bool adjoint = config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint(); + bool tkeNeeded = (rans && config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); int Unst_RestartIter = 0; unsigned long iPoint, iMarker, counter_local = 0, counter_global = 0; @@ -240,6 +241,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Density_Inf = config->GetDensity_FreeStreamND(); Energy_Inf = config->GetEnergy_FreeStreamND(); Mach_Inf = config->GetMach(); + const su2double TKE_Inf = config->GetTke_FreeStreamND(); /*--- Initialize the secondary values for direct derivative approxiations ---*/ @@ -307,6 +309,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Velocity2 += pow(nodes->GetSolution(iPoint,iDim+1)/Density,2); StaticEnergy= nodes->GetEnergy(iPoint) - 0.5*Velocity2; + if (tkeNeeded) StaticEnergy -= TKE_Inf; GetFluidModel()->SetTDState_rhoe(Density, StaticEnergy); Pressure= GetFluidModel()->GetPressure(); @@ -319,7 +322,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Solution[0] = Density_Inf; for (iDim = 0; iDim < nDim; iDim++) Solution[iDim+1] = Velocity_Inf[iDim]*Density_Inf; - Solution[nDim+1] = Energy_Inf*Density_Inf; + Solution[nDim+1] = Energy_Inf*Density_Inf; nodes->SetSolution(iPoint,Solution); nodes->SetSolution_Old(iPoint,Solution); counter_local++; @@ -801,7 +804,8 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes bool viscous = config->GetViscous(); bool gravity = config->GetGravityForce(); bool turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); - bool tkeNeeded = (turbulent && config->GetKind_Turb_Model() == TURB_MODEL::SST); + bool tkeNeeded = (turbulent && config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); + // bool tkeNeeded = (turbulent && config->GetKind_Turb_Model() == TURB_MODEL::SST); bool free_stream_temp = (config->GetKind_FreeStreamOption() == FREESTREAM_OPTION::TEMPERATURE_FS); bool reynolds_init = (config->GetKind_InitOption() == REYNOLDS); bool aeroelastic = config->GetAeroelastic_Simulation(); @@ -966,6 +970,19 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); + if (config->GetSSTParsedOptions().newBC) { + su2double Omega_Freestream = 10 * ModVel_FreeStream / config->GetLDomain(); + Tke_FreeStream = Omega_Freestream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + } else if (config->GetSSTParsedOptions().sust) { + su2double Omega_Freestream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + /*-- Compute the freestream energy. ---*/ + + if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; + } else { @@ -976,9 +993,8 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes } - /*-- Compute the freestream energy. ---*/ - - if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; config->SetEnergy_FreeStream(Energy_FreeStream); + /*-- Set the freestream energy. ---*/ + config->SetEnergy_FreeStream(Energy_FreeStream); /*--- Compute non dimensional quantities. By definition, Lref is one because we have converted the grid to meters. ---*/ @@ -4783,7 +4799,7 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; bool viscous = config->GetViscous(); - bool tkeNeeded = config->GetKind_Turb_Model() == TURB_MODEL::SST; + bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); CVariable* turbNodes = nullptr; if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); @@ -5029,7 +5045,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const bool viscous = config->GetViscous(), implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT), gravity = config->GetGravityForce(), - tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); CVariable* turbNodes = nullptr; if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); @@ -5077,7 +5093,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const auto Density_i = nodes->GetDensity(iPoint); const auto Energy_i = nodes->GetEnergy(iPoint); - const su2double StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + su2double StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + if (tkeNeeded) StaticEnergy_i -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -5158,7 +5175,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); + // if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0); break; case STATIC_SUPERSONIC_INFLOW_PD: @@ -5185,7 +5203,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); + // if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0); break; case DENSITY_VELOCITY: @@ -5218,6 +5237,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0); break; default: @@ -5422,7 +5442,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, /*--- Turbulent kinetic energy ---*/ - if (config->GetKind_Turb_Model() == TURB_MODEL::SST) + if (tkeNeeded) visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); @@ -7196,7 +7216,7 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, su2double Gas_Constant = config->GetGas_ConstantND(); string Marker_Tag = config->GetMarker_All_TagBound(val_marker); bool gravity = (config->GetGravityForce()); - bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); CVariable* turbNodes = nullptr; if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); @@ -7282,8 +7302,8 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = P_Exit/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); - // if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); + // if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Conservative variables, using the derived quantities ---*/ V_outlet[0] = Pressure / ( Gas_Constant * Density); @@ -7366,7 +7386,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con const su2double Gas_Constant = config->GetGas_ConstantND(); const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); - const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); CVariable* turbNodes = nullptr; if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index f2f1d6f4c4e..4156ae90a03 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -133,7 +133,7 @@ unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, cons unsigned long nonPhysicalPoints = 0; const TURB_MODEL turb_model = config->GetKind_Turb_Model(); - const bool tkeNeeded = (turb_model == TURB_MODEL::SST); + const bool tkeNeeded = (turb_model == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); AD::StartNoSharedReading(); @@ -799,8 +799,7 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, const unsigned short max_iter = config->GetwallModel_MaxIter(); const su2double relax = config->GetwallModel_RelFac(); - const bool SSTm = config->GetSSTParsedOptions().modified; - const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !SSTm; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); const auto* Node_Turb = (tkeNeeded) ? solver_container[TURB_SOL]->GetNodes() : nullptr; /*--- Compute the recovery factor From 317f5639781184a5ab19f794c955b5d3b687f535 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 20 Sep 2024 19:59:31 +0200 Subject: [PATCH 39/44] - Fixed tke integration in thermodynamic variables --- Common/include/CConfig.hpp | 2 +- .../numerics_simd/flow/convection/common.hpp | 7 +------ .../include/numerics_simd/flow/convection/roe.hpp | 8 ++------ SU2_CFD/include/numerics_simd/flow/variables.hpp | 2 +- SU2_CFD/src/solvers/CEulerSolver.cpp | 15 ++++++++++++--- SU2_CFD/src/solvers/CIncNSSolver.cpp | 2 +- 6 files changed, 18 insertions(+), 18 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index deac682f12f..7a46555c8a6 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1961,7 +1961,7 @@ class CConfig { * \brief Get the value of the non-dimensionalized freestream energy. * \return Non-dimensionalized freestream energy. */ - su2double GetEnergy_FreeStreamND(void) const { cout << "La chiedo non-dimensionale " << endl; return Energy_FreeStreamND; } + su2double GetEnergy_FreeStreamND(void) const { return Energy_FreeStreamND; } /*! * \brief Get the value of the non-dimensionalized freestream viscosity. diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index e3781a0c09d..0ca03b9a962 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -103,7 +103,6 @@ FORCEINLINE void musclEdgeLimited(Int iPoint, template FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int jPoint, bool muscl, LIMITER limiterType, - bool musclTurb, LIMITER limiterTypeTurb, const CPair& V1st, const VectorDbl& vector_ij, const VariableType& solution, @@ -113,9 +112,6 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int const auto& gradients = solution.GetGradient_Reconstruction(); const auto& limiters = solution.GetLimiter_Primitive(); - // const auto& gradientsTurb = turbSolution.GetGradient_Reconstruction(); - // const auto& limitersTurb = turbSolution.GetLimiter_Primitive(); - CPair V; for (size_t iVar = 0; iVar < ReconVarType::nVar; ++iVar) { @@ -153,10 +149,9 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int for (size_t iDim = 0; iDim < nDim; ++iDim) { v_squared += pow(R*V.j.velocity(iDim) + V.i.velocity(iDim), 2); } - Double tke = R*V1st.j.allTurb(0) + V1st.i.allTurb(0); /*--- Multiply enthalpy by R+1 since v^2 was not divided by (R+1)^2. * Note: a = sqrt((gamma-1) * (H - 0.5 * v^2)) ---*/ - const Double neg_sound_speed = enthalpy * (R+1) < (0.5 * v_squared - tke); + const Double neg_sound_speed = enthalpy * (R+1) < 0.5 * v_squared; /*--- Revert to first order if the state is non-physical. ---*/ Double bad_recon = fmax(neg_p_or_rho, neg_sound_speed); diff --git a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp index dabb24837cf..c0dbc2fb8ca 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp @@ -60,9 +60,7 @@ class CRoeBase : public Base { const bool finestGrid; const bool dynamicGrid; const bool muscl; - const bool musclTurb; const LIMITER typeLimiter; - const LIMITER typeLimiterTurb; using Base::turbVars; @@ -77,9 +75,7 @@ class CRoeBase : public Base { finestGrid(iMesh == MESH_0), dynamicGrid(config.GetDynamic_Grid()), muscl(finestGrid && config.GetMUSCL_Flow()), - musclTurb(finestGrid && config.GetMUSCL_Turb()), - typeLimiter(config.GetKind_SlopeLimit_Flow()), - typeLimiterTurb(config.GetKind_SlopeLimit_Turb()) { + typeLimiter(config.GetKind_SlopeLimit_Flow()) { } public: @@ -132,7 +128,7 @@ class CRoeBase : public Base { } auto V = reconstructPrimitives >( - iEdge, iPoint, jPoint, muscl, typeLimiter, musclTurb, typeLimiterTurb, V1st, vector_ij, solution, tkeNeeded); + iEdge, iPoint, jPoint, muscl, typeLimiter, V1st, vector_ij, solution, tkeNeeded); /*--- Compute conservative variables. ---*/ diff --git a/SU2_CFD/include/numerics_simd/flow/variables.hpp b/SU2_CFD/include/numerics_simd/flow/variables.hpp index 44fd815911d..028b26de2fe 100644 --- a/SU2_CFD/include/numerics_simd/flow/variables.hpp +++ b/SU2_CFD/include/numerics_simd/flow/variables.hpp @@ -136,7 +136,7 @@ FORCEINLINE CRoeVariables roeAveragedVariables(Double gamma, } roeAvg.enthalpy = (R*V.j.enthalpy() + V.i.enthalpy()) * D; roeAvg.tke = (R*V.j.tke() + V.i.tke()) * D; - roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity) - roeAvg.tke)); + roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity))); roeAvg.projVel = dot(roeAvg.velocity, normal); return roeAvg; } diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 1a2b494eea1..79a25cc18dc 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -5494,7 +5494,9 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain unsigned short nSpanWiseSections = geometry->GetnSpanWiseSections(config->GetMarker_All_TurbomachineryFlag(val_marker)); bool viscous = config->GetViscous(); bool gravity = (config->GetGravityForce()); - bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + bool tkeNeeded = ((config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified)); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); su2double *Normal, *turboNormal, *UnitNormal, *FlowDirMix, FlowDirMixMag, *turboVelocity; Normal = new su2double[nDim]; @@ -5568,6 +5570,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Energy_i = nodes->GetEnergy(iPoint); StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + if(tkeNeeded) StaticEnergy_i -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -5613,11 +5616,12 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain turboVelocity[iDim] = sqrt(Velocity2_e)*Flow_Dir[iDim]; ComputeBackVelocity(turboVelocity,turboNormal, Velocity_e, config->GetMarker_All_TurbomachineryFlag(val_marker),config->GetKind_TurboMachinery(iZone)); StaticEnthalpy_e = Enthalpy_e - 0.5 * Velocity2_e; + if(tkeNeeded) StaticEnthalpy_e -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_hs(StaticEnthalpy_e, Entropy_e); Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); + if(tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; case MIXING_IN: @@ -5648,10 +5652,12 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain ComputeBackVelocity(turboVelocity,turboNormal, Velocity_e, config->GetMarker_All_TurbomachineryFlag(val_marker),config->GetKind_TurboMachinery(iZone)); StaticEnthalpy_e = Enthalpy_e - 0.5 * Velocity2_e; + if(tkeNeeded) StaticEnthalpy_e -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_hs(StaticEnthalpy_e, Entropy_e); Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; + if(tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); // if (tkeNeeded) Energy_e += GetTke_Inf(); break; @@ -5670,6 +5676,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if(tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; case STATIC_PRESSURE: @@ -5687,6 +5694,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if(tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; @@ -5704,6 +5712,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if(tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; default: @@ -6915,7 +6924,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, const su2double Gas_Constant = config->GetGas_ConstantND(); const auto Kind_Inlet = config->GetKind_Inlet(); const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); - const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); /*--- Loop over all the vertices on this boundary marker ---*/ diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index d8327f1297a..e5a24427eaa 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -287,7 +287,7 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, c const TURB_MODEL turb_model = config->GetKind_Turb_Model(); const SPECIES_MODEL species_model = config->GetKind_Species_Model(); - bool tkeNeeded = (turb_model == TURB_MODEL::SST); + bool tkeNeeded = (turb_model == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); AD::StartNoSharedReading(); From d6af7021814671fdcda38c569d0ba08079136e15 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 2 Oct 2024 19:15:45 +0200 Subject: [PATCH 40/44] - Code improvements --- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 8c3984b2f51..a8a7788353b 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -557,7 +557,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; - const bool SSTm = config->GetSSTParsedOptions().modified; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified); /*--- Normalized normal vector ---*/ @@ -593,7 +593,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); Mean_Thermal_Conductivity = 0.5*(Thermal_Conductivity_i + Thermal_Conductivity_j); - const su2double Mean_turb_ke_ForStressTensor = !SSTm ? Mean_turb_ke : 0.0; + const su2double Mean_turb_ke_ForStressTensor = tkeNeeded ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ @@ -877,7 +877,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; - const bool SSTm = config->GetSSTParsedOptions().modified; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified); /*--- Normalized normal vector ---*/ @@ -923,7 +923,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); Mean_Thermal_Conductivity = 0.5*(Thermal_Conductivity_i + Thermal_Conductivity_j); Mean_Cp = 0.5*(Cp_i + Cp_j); - const su2double Mean_turb_ke_ForStressTensor = !SSTm ? Mean_turb_ke : 0.0; + const su2double Mean_turb_ke_ForStressTensor = tkeNeeded ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ From 8f97a162c09ba02d6bcd4e1bcbca07cd0b06e472 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 11 Oct 2024 11:56:09 +0200 Subject: [PATCH 41/44] - minor changes --- Common/include/option_structure.hpp | 1 - SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 2 +- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index b8be96692b1..e7cab6e7bd3 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1002,7 +1002,6 @@ static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) MakePair("V1994m", SST_OPTIONS::V1994m) MakePair("V2003m", SST_OPTIONS::V2003m) - /// TODO: For now we do not support "unmodified" versions of SST. MakePair("V1994", SST_OPTIONS::V1994) MakePair("V2003", SST_OPTIONS::V2003) MakePair("SUSTAINING", SST_OPTIONS::SUST) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index be2615c5e2f..c1478f2e2ab 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -880,7 +880,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { if ( P > prod_limit ) PLim = 1.0; su2double pk = max(0.0, min(P, prod_limit)); - const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; + const su2double eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; const su2double zeta = max(ScalarVar_i[1], eddy_visc_var * F2_i / a1); /*--- Production limiter only for V2003, recompute for V1994. ---*/ diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 58fbee558b1..59f2e2593b6 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -251,7 +251,7 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai const su2double kine = nodes->GetSolution(iPoint,0); const su2double omega = nodes->GetSolution(iPoint,1); - const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag; + const su2double eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag; const su2double muT = max(0.0, rho * a1 * kine / max(a1 * omega, eddy_visc_var * F2)); nodes->SetmuT(iPoint, muT); From a8a603e37a431bf6c09fb35139b3bf1f6abab959 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 11 Oct 2024 12:39:10 +0200 Subject: [PATCH 42/44] - update externals --- SU2_CFD/include/numerics_simd/flow/convection/roe.hpp | 2 -- externals/codi | 2 +- externals/medi | 2 +- externals/opdi | 2 +- 4 files changed, 3 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp index fdf11e7bec4..98be310304c 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp @@ -105,8 +105,6 @@ class CRoeBase : public Base { const auto iPoint = geometry.edges->GetNode(iEdge,0); const auto jPoint = geometry.edges->GetNode(iEdge,1); - // cout << gatherVariables(iPoint, turbVars->GetSolution()) << endl; - /*--- Geometric properties. ---*/ const auto vector_ij = distanceVector(iPoint, jPoint, geometry.nodes->GetCoord()); diff --git a/externals/codi b/externals/codi index 762ba7698e3..c6b039e5c9e 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit 762ba7698e3ceaa1b17aa299421ddd418f00b823 +Subproject commit c6b039e5c9edb7675f90ffc725f9dd8e66571264 diff --git a/externals/medi b/externals/medi index 7d550831e0e..ab3a7688f6d 160000 --- a/externals/medi +++ b/externals/medi @@ -1 +1 @@ -Subproject commit 7d550831e0e233a85b9d9af9c181d7ecb2929946 +Subproject commit ab3a7688f6d518f8d940eb61a341d89f51922ba4 diff --git a/externals/opdi b/externals/opdi index a6b9655c240..8c897988172 160000 --- a/externals/opdi +++ b/externals/opdi @@ -1 +1 @@ -Subproject commit a6b9655c240af2a35454a61727e5bbbbaa3a425f +Subproject commit 8c89798817253abb017d857a0ae7f0520187645c From 340faa48533856fdcba6b897de5dd42ae4df61e4 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 11 Oct 2024 18:55:46 +0200 Subject: [PATCH 43/44] - Fixed inlet BCs when new BCs are used --- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 991830a5f3b..def030f630c 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -645,8 +645,7 @@ void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C const su2double VelMag2 = GeometryToolbox::SquaredNorm(nDim, Velocity_Inlet); if (sstParsedOptions.newBC) { - su2double LDomain = 0.53; - Inlet_Vars[1] = 10 * sqrt(VelMag2) / LDomain; + Inlet_Vars[1] = 10 * sqrt(VelMag2) / config->GetLDomain(); Inlet_Vars[0] = Inlet_Vars[1]*(Laminar_Viscosity_Inlet*viscRatio)/Density_Inlet; } else { Inlet_Vars[0] = 3.0 / 2.0 * (VelMag2 * pow(Intensity, 2)); From c4f10fef4fbec3012d8cdb7f5813ef210f1a8cff Mon Sep 17 00:00:00 2001 From: rois1995 <43813346+rois1995@users.noreply.github.com> Date: Thu, 31 Oct 2024 09:53:29 +0100 Subject: [PATCH 44/44] Update Common/include/CConfig.hpp Co-authored-by: Nijso --- Common/include/CConfig.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index abb375fa9a4..ee47d4cacaf 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -865,7 +865,7 @@ class CConfig { nMu_Temperature_Ref, /*!< \brief Number of species reference temperature for Sutherland model. */ nMu_S, /*!< \brief Number of species reference S for Sutherland model. */ nThermal_Conductivity_Constant,/*!< \brief Number of species constant thermal conductivity. */ - nPrandtl_Lam, /*!< \brief Prandtl number. */ + nPrandtl_Lam, /*!< \brief Number of species laminar Prandtl number. */ nPrandtl_Turb, /*!< \brief Number of species turbulent Prandtl number. */ nConstant_Lewis_Number; /*!< \brief Number of species Lewis Number. */ su2double Diffusivity_Constant; /*!< \brief Constant mass diffusivity for scalar transport. */