From 388f85a3da092d3b74be4786b2c9bdf692ad4151 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 15 Sep 2023 09:21:47 +0200 Subject: [PATCH 01/16] - Added computation of desired grid size for Scale Resolving Simulations based on k and w - Added computation of k and w starting from SA variables - Added first implementation of SAS (Scale Adaptive Simulation) model --- Common/include/option_structure.hpp | 5 +++ SU2_CFD/include/numerics/CNumerics.hpp | 7 ++++ .../numerics/turbulent/turb_sources.hpp | 20 +++++++-- SU2_CFD/include/variables/CTurbSAVariable.hpp | 24 +++++++++++ .../include/variables/CTurbSSTVariable.hpp | 14 +++++++ SU2_CFD/include/variables/CTurbVariable.hpp | 19 ++++++++- SU2_CFD/include/variables/CVariable.hpp | 41 +++++++++++++++++++ SU2_CFD/src/output/CFlowOutput.cpp | 6 +++ SU2_CFD/src/solvers/CTurbSASolver.cpp | 29 ++++++++++++- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 11 +++++ SU2_CFD/src/variables/CTurbSSTVariable.cpp | 2 + SU2_CFD/src/variables/CTurbVariable.cpp | 2 + 12 files changed, 173 insertions(+), 7 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index c8181a3a13f..dabddf8a974 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -990,6 +990,7 @@ 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. */ + SAS, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) @@ -1002,6 +1003,7 @@ static const MapType SST_Options_Map = { MakePair("VORTICITY", SST_OPTIONS::V) MakePair("KATO-LAUNDER", SST_OPTIONS::KL) MakePair("UQ", SST_OPTIONS::UQ) + MakePair("SAS", SST_OPTIONS::SAS) }; /*! @@ -1012,6 +1014,7 @@ struct SST_ParsedOptions { SST_OPTIONS production = SST_OPTIONS::NONE; /*!< \brief Enum for production corrections/modifiers for SST model. */ bool sust = false; /*!< \brief Bool for SST model with sustaining terms. */ bool uq = false; /*!< \brief Bool for using uncertainty quantification. */ + bool sas = false; /*!< \brief Bool for using Scale Adaptive Simulations. */ bool modified = false; /*!< \brief Bool for modified (m) SST model. */ }; @@ -1048,6 +1051,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne const bool sst_v = IsPresent(SST_OPTIONS::V); const bool sst_kl = IsPresent(SST_OPTIONS::KL); const bool sst_uq = IsPresent(SST_OPTIONS::UQ); + const bool sst_sas = IsPresent(SST_OPTIONS::SAS); if (sst_1994 && sst_2003) { SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION); @@ -1071,6 +1075,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.sust = sst_sust; SSTParsedOptions.modified = sst_m; SSTParsedOptions.uq = sst_uq; + SSTParsedOptions.sas = sst_sas; return SSTParsedOptions; } diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 14088f021cd..e90dae26a23 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 FTrans; /*!< \brief SAS function */ + public: /*! * \brief Return type used in some "ComputeResidual" overloads to give a @@ -707,6 +709,11 @@ class CNumerics { */ virtual void SetCrossDiff(su2double val_CDkw_i) {/* empty */}; + /*! + * \brief Get the value of the value of FTrans. + */ + inline virtual su2double GetFTrans() const { return 0.0; } + /*! * \brief Set the value of the effective intermittency for the LM model. * \param[in] intermittency_eff_i - Value of the effective intermittency at point i. diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index af920cd34c1..1b7d3baed1e 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -605,6 +605,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Closure constants ---*/ const su2double sigma_k_1, sigma_k_2, sigma_w_1, sigma_w_2, beta_1, beta_2, beta_star, a1, alfa_1, alfa_2; const su2double prod_lim_const; + const su2double cTrans; /*--- Ambient values for SST-SUST. ---*/ const su2double kAmb, omegaAmb; @@ -689,6 +690,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { alfa_1(constants[8]), alfa_2(constants[9]), prod_lim_const(constants[10]), + cTrans(1.25), kAmb(val_kine_Inf), omegaAmb(val_omega_Inf) { /*--- "Allocate" the Jacobian using the static buffer. ---*/ @@ -848,8 +850,11 @@ class CSourcePieceWise_TurbSST final : public CNumerics { Residual[1] += pw * Volume; /*--- Add the dissipation terms to the residuals.---*/ - - Residual[0] -= dk * Volume; + FTrans = 1.0; + if (sstParsedOptions.sas) { + FTrans = max(1.0, pow(StrainMag_i / (cTrans * VorticityMag), 2.0)); + } + Residual[0] -= dk * Volume * FTrans; Residual[1] -= dw * Volume; /*--- Cross diffusion ---*/ @@ -862,8 +867,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Implicit part ---*/ - Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume; - Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume; + Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * FTrans; + Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * FTrans; Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume; } @@ -873,4 +878,11 @@ class CSourcePieceWise_TurbSST final : public CNumerics { return ResidualType<>(Residual, Jacobian_i, nullptr); } + + /*! + * \brief Get the value of the FTrans. + */ + inline su2double GetFTrans() const override { return FTrans; } + + }; diff --git a/SU2_CFD/include/variables/CTurbSAVariable.hpp b/SU2_CFD/include/variables/CTurbSAVariable.hpp index 569e416fa15..b76539d4c05 100644 --- a/SU2_CFD/include/variables/CTurbSAVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp @@ -42,6 +42,8 @@ class CTurbSAVariable final : public CTurbVariable { VectorType DES_LengthScale; VectorType Vortex_Tilting; + VectorType k, Omega; /*!< \brief SST variables as computed through SA solution. */ + public: /*! * \brief Constructor of the class. @@ -87,4 +89,26 @@ class CTurbSAVariable final : public CTurbVariable { */ inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); } + /*! + * \brief Get the value of the turbulence kinetic energy. + * \return the value of the turbulence kinetic energy. + */ + inline su2double GetSSTVariables_k(unsigned long iPoint) const { return k(iPoint); } + + /*! + * \brief Get the value of the turbulence frequency Omega. + * \return the value of the turbulence frequency Omega. + */ + inline su2double GetSSTVariables_omega(unsigned long iPoint) const { return Omega(iPoint); } + + /*! + * \brief Set the value of the SST variables computed with SA solution. + * \param[in] val_k + * \param[in] val_Omega + */ + void SetSSTVariables(unsigned long iPoint, su2double val_k, su2double val_Omega) { + k(iPoint) = val_k; + Omega(iPoint) = val_Omega; + } + }; diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp index 0113a0ead3a..384fe86f95a 100644 --- a/SU2_CFD/include/variables/CTurbSSTVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSSTVariable.hpp @@ -45,6 +45,7 @@ class CTurbSSTVariable final : public CTurbVariable { VectorType F2; /*!< \brief Menter blending function for blending of k-w and k-eps. */ VectorType CDkw; /*!< \brief Cross-diffusion. */ SST_ParsedOptions sstParsedOptions; + VectorType FTrans; public: /*! * \brief Constructor of the class. @@ -87,4 +88,17 @@ class CTurbSSTVariable final : public CTurbVariable { * \brief Get the value of the cross diffusion of tke and omega. */ inline su2double GetCrossDiff(unsigned long iPoint) const override { return CDkw(iPoint); } + + /*! + * \brief Get the value of the FTrans. + */ + inline su2double GetFTrans(unsigned long iPoint) const override { return FTrans(iPoint); } + + /*! + * \brief Set the value of the FTrans. + */ + inline void SetFTrans(unsigned long iPoint, su2double val_FTrans) override { + FTrans(iPoint) = val_FTrans; + }; + }; diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index 1a476fdac70..ba3a0651253 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -41,8 +41,9 @@ class CTurbVariable : public CScalarVariable { public: static constexpr size_t MAXNVAR = 2; - VectorType turb_index; - VectorType intermittency; /*!< \brief Value of the intermittency for the trans. model. */ + VectorType turb_index; /*!< \brief Value of the turbulence index for transition simulations. */ + VectorType intermittency; /*!< \brief Value of the intermittency for the transition model. */ + VectorType SRSGridSize; /*!< \brief alue of the desired grid size for Scale Resolving Simulations. */ /*! * \brief Constructor of the class. @@ -99,6 +100,20 @@ class CTurbVariable : public CScalarVariable { * \param[in] val_intermittency - New value of the intermittency. */ inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; } + + /*! + * \brief Get the desired grid size for Scale Resolving Simulations. + * \param[in] iPoint - Point index. + * \return the value of the desired grid size for Scale Resolving Simulations. + */ + inline su2double GetSRSGridSize(unsigned long iPoint) const final { return SRSGridSize(iPoint); } + + /*! + * \brief Set the value of the desired grid size for Scale Resolving Simulations. + * \param[in] iPoint - Point index. + * \param[in] val_muT - Value of the desired grid size for Scale Resolving Simulations. + */ + inline void SetSRSGridSize(unsigned long iPoint, su2double val_gridSize) final { SRSGridSize(iPoint) = val_gridSize; } }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 050be2b7787..363fad7a7b7 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1647,6 +1647,16 @@ class CVariable { */ inline virtual su2double GetCrossDiff(unsigned long iPoint) const { return 0.0; } + /*! + * \brief Get the value of the value of FTrans. + */ + inline virtual su2double GetFTrans(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief Set the value of the value of FTrans. + */ + inline virtual void SetFTrans(unsigned long iPoint, su2double val_FTrans) {} + /*! * \brief Get the value of the eddy viscosity. * \return the value of the eddy viscosity. @@ -1695,6 +1705,37 @@ class CVariable { */ inline virtual void SetmuT(unsigned long iPoint, su2double val_muT) {} + /*! + * \brief Get the value of the turbulence kinetic energy. + * \return the value of the turbulence kinetic energy. + */ + inline virtual su2double GetSSTVariables_k(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief Get the value of the Rate of dissipation Omega. + * \return the value of the Rate of dissipation Omega. + */ + inline virtual su2double GetSSTVariables_omega(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief Set the value of the SST variables computed with SA solution. + * \param[in] val_k + * \param[in] val_Omega + */ + inline virtual void SetSSTVariables(unsigned long iPoint, su2double val_k, su2double val_Omega) {} + + /*! + * \brief Get the value of the desired grid size for Scale Resolving Simulations. + * \return the value of the desired grid size for Scale Resolving Simulations. + */ + inline virtual su2double GetSRSGridSize(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief Set the value of the desired grid size for Scale Resolving Simulations. + * \param[in] val_gridSize + */ + inline virtual void SetSRSGridSize(unsigned long iPoint, su2double val_gridSize) {} + /*! * \brief Set the value of the turbulence index. * \param[in] val_turb_index - turbulence index diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 75bd80feebb..bf9d989d202 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1243,11 +1243,14 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ switch (TurbModelFamily(config->GetKind_Turb_Model())) { case TURB_FAMILY::SA: AddVolumeOutput("NU_TILDE", "Nu_Tilde", "SOLUTION", "Spalart-Allmaras variable"); + AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations"); break; case TURB_FAMILY::KW: AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy"); AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); + AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations"); + if (config->GetSSTParsedOptions().sas) AddVolumeOutput("FTRANS", "FTrans", "SOLUTION", "value of FTrans for SAS simulation"); break; case TURB_FAMILY::NONE: @@ -1516,6 +1519,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con case TURB_FAMILY::SA: SetVolumeOutputValue("NU_TILDE", iPoint, Node_Turb->GetSolution(iPoint, 0)); SetVolumeOutputValue("RES_NU_TILDE", iPoint, turb_solver->LinSysRes(iPoint, 0)); + SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, Node_Turb->GetSRSGridSize(iPoint)); if (limiter) { SetVolumeOutputValue("LIMITER_NU_TILDE", iPoint, Node_Turb->GetLimiter(iPoint, 0)); } @@ -1526,6 +1530,8 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("DISSIPATION", iPoint, Node_Turb->GetSolution(iPoint, 1)); SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); + SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, Node_Turb->GetSRSGridSize(iPoint)); + if (config->GetSSTParsedOptions().sas) SetVolumeOutputValue("FTRANS", iPoint, Node_Turb->GetFTrans(iPoint)); if (limiter) { SetVolumeOutputValue("LIMITER_TKE", iPoint, Node_Turb->GetLimiter(iPoint, 0)); SetVolumeOutputValue("LIMITER_DISSIPATION", iPoint, Node_Turb->GetLimiter(iPoint, 1)); diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 2d301dc7722..a81f75a5e09 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -224,7 +224,8 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain const su2double nu = mu/rho; const su2double nu_hat = nodes->GetSolution(iPoint,0); const su2double roughness = geometry->nodes->GetRoughnessHeight(iPoint); - const su2double dist = geometry->nodes->GetWall_Distance(iPoint) + rough_const * roughness; + const su2double wallDistance = geometry->nodes->GetWall_Distance(iPoint); + const su2double dist = wallDistance + rough_const * roughness; su2double Ji = nu_hat/nu; if (roughness > 1.0e-10) @@ -239,6 +240,32 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain nodes->SetmuT(iPoint,muT); + // Now compute k and w + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + + const su2double BetaStar = 0.09; + const su2double a1 = 0.31; + + const su2double Omega = GeometryToolbox::Norm(nDim, flowNodes->GetVorticity(iPoint))/sqrt(BetaStar); + + const su2double nut = muT / rho; + const su2double arg2_1 = 2.0 * sqrt(nut / Omega) / (BetaStar * wallDistance); + const su2double arg2_2 = 500.0 * flowNodes->GetLaminarViscosity(iPoint) / (wallDistance*wallDistance * Omega); + + const su2double arg2 = max(arg2_1, arg2_2); + const su2double F2 = tanh(arg2*arg2); + + const su2double k = nut * max(Omega, flowNodes->GetStrainMag(iPoint) * F2 / a1); + + nodes->SetSSTVariables(iPoint, k, Omega); + + // Now compute desired cell size for Scale Resolving Simulations + const su2double Cmu = 0.09; + const su2double RANSLength = sqrt(k) / max(1e-20, (Cmu * Omega)); + const su2double RatioL = 0.1; // it should be less or equal than 0.2 - 0.1. Should be taken as input from config? + const su2double SRSGridSize = RANSLength * RatioL; + nodes->SetSRSGridSize(iPoint, SRSGridSize); + } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index e1e9ef644c7..e149df25bd8 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -236,6 +236,12 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai nodes->SetmuT(iPoint, muT); + // Now compute desired cell size for Scale Resolving Simulations + const su2double RANSLength = sqrt(nodes->GetSolution(iPoint, 0)) / max(1e-20, (constants[6] * nodes->GetSolution(iPoint, 1))); + const su2double RatioL = 0.1; // it should be less or equal than 0.2 - 0.1. Should be taken as input from config? + const su2double SRSGridSize = RANSLength * RatioL; + nodes->SetSRSGridSize(iPoint, SRSGridSize); + } END_SU2_OMP_FOR @@ -362,6 +368,11 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta auto residual = numerics->ComputeResidual(config); + /*--- Store the SAS function ---*/ + if (sstParsedOptions.sas) { + nodes->SetFTrans(iPoint, numerics->GetFTrans()); + } + /*--- Store the intermittency ---*/ if (config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) { diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index 9778f1e37cc..4a477665f23 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -51,6 +51,8 @@ CTurbSSTVariable::CTurbSSTVariable(su2double kine, su2double omega, su2double mu CDkw.resize(nPoint) = su2double(0.0); muT.resize(nPoint) = mut; + + if(sstParsedOptions.sas) FTrans.resize(nPoint) = su2double(1.0); } void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_viscosity, diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp index 54bede5033f..2b2069b3590 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -35,4 +35,6 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned turb_index.resize(nPoint) = su2double(1.0); intermittency.resize(nPoint) = su2double(1.0); + SRSGridSize.resize(nPoint) = su2double(0.0); + } From 41cb27e83c1d19b475644f82a639321dacbd27d7 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 15 Sep 2023 16:25:24 +0200 Subject: [PATCH 02/16] - Added a more complex SAS model named SAS_COMPLICATED - Changed previous SAS model to SAS_SIMPLE --- Common/include/option_structure.hpp | 20 +++++- SU2_CFD/include/numerics/CNumerics.hpp | 16 +++++ .../numerics/turbulent/turb_sources.hpp | 47 +++++++++++- .../include/variables/CTurbSSTVariable.hpp | 49 +++++++++++++ SU2_CFD/include/variables/CVariable.hpp | 30 ++++++++ SU2_CFD/src/solvers/CSolver.cpp | 71 +++++++++++++++++++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 62 +++++++++++++++- SU2_CFD/src/variables/CTurbSSTVariable.cpp | 7 +- 8 files changed, 296 insertions(+), 6 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index dabddf8a974..e00ac40d7dd 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -990,7 +990,8 @@ 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. */ - SAS, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ + SAS_SIMPLE, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ + SAS_COMPLICATED, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) @@ -1003,7 +1004,8 @@ static const MapType SST_Options_Map = { MakePair("VORTICITY", SST_OPTIONS::V) MakePair("KATO-LAUNDER", SST_OPTIONS::KL) MakePair("UQ", SST_OPTIONS::UQ) - MakePair("SAS", SST_OPTIONS::SAS) + MakePair("SAS_SIMPLE", SST_OPTIONS::SAS_SIMPLE) + MakePair("SAS_COMPLICATED", SST_OPTIONS::SAS_COMPLICATED) }; /*! @@ -1012,6 +1014,7 @@ static const MapType SST_Options_Map = { struct SST_ParsedOptions { SST_OPTIONS version = SST_OPTIONS::V1994; /*!< \brief Enum SST base model. */ SST_OPTIONS production = SST_OPTIONS::NONE; /*!< \brief Enum for production corrections/modifiers for SST model. */ + SST_OPTIONS sasModel = SST_OPTIONS::SAS_SIMPLE; /*!< \brief Enum SST base model. */ bool sust = false; /*!< \brief Bool for SST model with sustaining terms. */ bool uq = false; /*!< \brief Bool for using uncertainty quantification. */ bool sas = false; /*!< \brief Bool for using Scale Adaptive Simulations. */ @@ -1051,7 +1054,16 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne const bool sst_v = IsPresent(SST_OPTIONS::V); const bool sst_kl = IsPresent(SST_OPTIONS::KL); const bool sst_uq = IsPresent(SST_OPTIONS::UQ); - const bool sst_sas = IsPresent(SST_OPTIONS::SAS); + const bool sst_sas_simple = IsPresent(SST_OPTIONS::SAS_SIMPLE); + const bool sst_sas_comp = IsPresent(SST_OPTIONS::SAS_COMPLICATED); + const bool sst_sas = sst_sas_simple || sst_sas_comp; + if (sst_sas_simple && sst_sas_comp) { + SU2_MPI::Error("Two versions (Simple and Complicated) selected for SAS under SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION); + } else if (sst_sas_simple) { + SSTParsedOptions.sasModel = SST_OPTIONS::SAS_SIMPLE; + } else { + SSTParsedOptions.sasModel = SST_OPTIONS::SAS_COMPLICATED; + } if (sst_1994 && sst_2003) { SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION); @@ -2455,6 +2467,7 @@ enum PERIODIC_QUANTITIES { PERIODIC_LIM_PRIM_1 , /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */ PERIODIC_LIM_PRIM_2 , /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */ PERIODIC_IMPLICIT , /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */ + PERIODIC_VEL_LAPLACIAN , /*!< \brief Velocity Laplacian communication for SAS (periodic only). */ }; /*! @@ -2486,6 +2499,7 @@ enum MPI_QUANTITIES { MESH_DISPLACEMENTS , /*!< \brief Mesh displacements at the interface. */ SOLUTION_TIME_N , /*!< \brief Solution at time n. */ SOLUTION_TIME_N1 , /*!< \brief Solution at time n-1. */ + VELOCITY_LAPLACIAN , /*!< \brief Velocity Laplacian communication. */ }; /*! diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index e90dae26a23..5b375d05887 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -191,6 +191,7 @@ class CNumerics { bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ su2double FTrans; /*!< \brief SAS function */ + su2double VelLapl_X, VelLapl_Y, VelLapl_Z; public: /*! @@ -714,6 +715,21 @@ class CNumerics { */ inline virtual su2double GetFTrans() const { return 0.0; } + /*! + * \brief Get the value of the value of FTrans. + */ + inline void SetVelLapl(su2double val_VelLapl_X, su2double val_VelLapl_Y) { + VelLapl_X = val_VelLapl_X; + VelLapl_Y = val_VelLapl_Y; + } + + /*! + * \brief Get the value of the value of FTrans. + */ + inline void SetVelLapl_Z(su2double val_VelLapl_Z) { + VelLapl_Z = val_VelLapl_Z; + } + /*! * \brief Set the value of the effective intermittency for the LM model. * \param[in] intermittency_eff_i - Value of the effective intermittency at point i. diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 1b7d3baed1e..87e6657e053 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -849,9 +849,50 @@ class CSourcePieceWise_TurbSST final : public CNumerics { Residual[0] += pk * Volume; Residual[1] += pw * Volume; + /*--- Compute SAS source term. ---*/ + su2double Q_SAS = 0.0; + + if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED) { + + const su2double KolmConst = 0.41; + const su2double csi2 = 3.51; + const su2double sigma_phi = 2.0/3.0; + const su2double C = 2.0; + const su2double C_s = 0.5; // Honestly I do not know if it is the right one. + const su2double gridSize = pow(Volume, 1.0/3.0); + // Scale of the modeled turbulence + const su2double L = sqrt(ScalarVar_i[0]) / (pow(beta_star, 0.25) * ScalarVar_i[1]); + // Von Karman Length Scale + su2double VelLaplMag = VelLapl_X*VelLapl_X + VelLapl_Y*VelLapl_Y; + if (nDim == 3) VelLaplMag += VelLapl_Z*VelLapl_Z; + const su2double L_vK_1 = KolmConst * StrainMag_i / sqrt(VelLaplMag); + const su2double L_vK_2 = C_s * sqrt(KolmConst * csi2 / (beta_blended/beta_star - alfa_blended)) * gridSize; + const su2double L_vK = max(L_vK_1, L_vK_2); + + su2double gradOmega = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + gradOmega += ScalarVar_Grad_i[1][iDim]*ScalarVar_Grad_i[1][iDim]; + + gradOmega /= ScalarVar_i[1]; + + su2double gradTKE = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + gradTKE += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim]; + + gradTKE /= ScalarVar_i[0]; + + const su2double Q_SAS_1 = csi2 * KolmConst * StrainMag_i * StrainMag_i * (L/L_vK) * (L/L_vK); + const su2double Q_SAS_2 = C * (2*ScalarVar_i[0] / sigma_phi) * max(gradOmega, gradTKE); + + Q_SAS = max(Q_SAS_1 - Q_SAS_2, 0.0); + + Residual[1] += Density_i * Q_SAS * Volume; + + } + /*--- Add the dissipation terms to the residuals.---*/ FTrans = 1.0; - if (sstParsedOptions.sas) { + if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_SIMPLE) { FTrans = max(1.0, pow(StrainMag_i / (cTrans * VorticityMag), 2.0)); } Residual[0] -= dk * Volume * FTrans; @@ -871,6 +912,10 @@ class CSourcePieceWise_TurbSST final : public CNumerics { Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * FTrans; Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume; + + if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED) { + Jacobian_i[0][0] += Q_SAS * Volume / ScalarVar_i[0]; + } } AD::SetPreaccOut(Residual, nVar); diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp index 384fe86f95a..5cfe6270569 100644 --- a/SU2_CFD/include/variables/CTurbSSTVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSSTVariable.hpp @@ -46,6 +46,10 @@ class CTurbSSTVariable final : public CTurbVariable { VectorType CDkw; /*!< \brief Cross-diffusion. */ SST_ParsedOptions sstParsedOptions; VectorType FTrans; + VectorType VelocityLaplacian_X; + VectorType VelocityLaplacian_Y; + VectorType VelocityLaplacian_Z; + public: /*! * \brief Constructor of the class. @@ -101,4 +105,49 @@ class CTurbSSTVariable final : public CTurbVariable { FTrans(iPoint) = val_FTrans; }; + /*! + * \brief Get the value of the FTrans. + */ + inline su2double GetVelLapl_X(unsigned long iPoint) const override { return VelocityLaplacian_X(iPoint); } + + /*! + * \brief Get the value of the FTrans. + */ + inline su2double GetVelLapl_Y(unsigned long iPoint) const override { return VelocityLaplacian_Y(iPoint); } + + /*! + * \brief Get the value of the FTrans. + */ + inline su2double GetVelLapl_Z(unsigned long iPoint) const override { return VelocityLaplacian_Z(iPoint); } + + /*! + * \brief Set the value of the FTrans. + */ + inline void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) override { + VelocityLaplacian_X(iPoint) += val_VelLapl_X; + VelocityLaplacian_Y(iPoint) += val_VelLapl_Y; + }; + + /*! + * \brief Set the value of the FTrans. + */ + inline void AddVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) override { + VelocityLaplacian_Z(iPoint) += val_VelLapl_Z; + }; + + /*! + * \brief Set the value of the FTrans. + */ + inline void SetVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) override { + VelocityLaplacian_X(iPoint) = val_VelLapl_X; + VelocityLaplacian_Y(iPoint) = val_VelLapl_Y; + }; + + /*! + * \brief Set the value of the FTrans. + */ + inline void SetVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) override { + VelocityLaplacian_Z(iPoint) = val_VelLapl_Z; + }; + }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 363fad7a7b7..82fe43caa47 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1657,6 +1657,36 @@ class CVariable { */ inline virtual void SetFTrans(unsigned long iPoint, su2double val_FTrans) {} + /*! + * \brief Get the value of the value of FTrans. + */ + inline virtual su2double GetVelLapl_X(unsigned long iPoint) const { return 0.0; } + /*! + * \brief Get the value of the value of FTrans. + */ + inline virtual su2double GetVelLapl_Y(unsigned long iPoint) const { return 0.0; } + /*! + * \brief Get the value of the value of FTrans. + */ + inline virtual su2double GetVelLapl_Z(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief Set the value of the value of FTrans. + */ + inline virtual void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) {} + /*! + * \brief Set the value of the value of FTrans. + */ + inline virtual void AddVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) {} + /*! + * \brief Set the value of the value of FTrans. + */ + inline virtual void SetVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) {} + /*! + * \brief Set the value of the value of FTrans. + */ + inline virtual void SetVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) {} + /*! * \brief Get the value of the eddy viscosity. * \return the value of the eddy viscosity. diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 5cd2fd93425..a0afe86738e 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -364,6 +364,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, auto *Diff = new su2double[nVar]; auto *Und_Lapl = new su2double[nVar]; + auto *Vel_Lapl = new su2double[nDim]; auto *Sol_Min = new su2double[nPrimVarGrad]; auto *Sol_Max = new su2double[nPrimVarGrad]; auto *rotPrim_i = new su2double[nPrimVar]; @@ -969,6 +970,58 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, } break; + case PERIODIC_VEL_LAPLACIAN: + + /*--- For JST, the undivided Laplacian must be computed + consistently by using the complete control volume info + from both sides of the periodic face. ---*/ + + for (iDim = 0; iDim < nDim; iDim++) + Vel_Lapl[iDim] = 0.0; + + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + + /*--- Avoid periodic boundary points so that we do not + duplicate edges on both sides of the periodic BC. ---*/ + + if (!geometry->nodes->GetPeriodicBoundary(jPoint)) { + + /*--- Solution differences ---*/ + const su2double distance = GeometryToolbox::Distance(nDim, base_nodes->GetMesh_Coord(iPoint), base_nodes->GetMesh_Coord(jPoint)); + + for (iDim = 0; iDim < nDim; iDim++) + Diff[iDim] = (base_nodes->GetVelocity(iPoint, iDim) - + base_nodes->GetVelocity(jPoint, iDim))/distance; + + + boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint); + boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint); + + /*--- Both points inside the domain, or both in the boundary ---*/ + /*--- iPoint inside the domain, jPoint on the boundary ---*/ + + if (!boundary_i || boundary_j) { + if (geometry->nodes->GetDomain(iPoint)){ + for (iDim = 0; iDim< nDim; iDim++) + Vel_Lapl[iDim] -= Diff[iDim]; + } + } + } + } + + /*--- Store the components to be communicated in the buffer. ---*/ + + for (iDim = 0; iDim < nDim; iDim++) + bufDSend[buf_offset+iDim] = Vel_Lapl[iDim]; + + /*--- Rotate the momentum components of the Laplacian. ---*/ + + if (rotate_periodic) { + Rotate(zeros, &Vel_Lapl[0], &bufDSend[buf_offset+0]); + } + + break; + default: SU2_MPI::Error("Unrecognized quantity for periodic communication.", @@ -1291,6 +1344,15 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, limiter(iPoint, iVar) = min(limiter(iPoint, iVar), bufDRecv[buf_offset+iVar]); break; + case PERIODIC_VEL_LAPLACIAN: + + /*--- Adjust the undivided Laplacian. The accumulation was + with a subtraction before communicating, so now just add. ---*/ + + base_nodes->AddVelLapl(iPoint, bufDRecv[buf_offset+0], bufDRecv[buf_offset+1]); + if(nDim == 3) base_nodes->AddVelLapl_Z(iPoint, bufDRecv[buf_offset+2]); + + break; default: @@ -1529,6 +1591,11 @@ void CSolver::InitiateComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution_time_n1(iPoint, iVar); break; + case VELOCITY_LAPLACIAN: + bufDSend[buf_offset+0] = base_nodes->GetVelLapl_X(iPoint); + bufDSend[buf_offset+1] = base_nodes->GetVelLapl_Y(iPoint); + if(nDim == 3) bufDSend[buf_offset+2] = base_nodes->GetVelLapl_Z(iPoint); + break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION); @@ -1677,6 +1744,10 @@ void CSolver::CompleteComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) base_nodes->Set_Solution_time_n1(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; + case VELOCITY_LAPLACIAN: + base_nodes->SetVelLapl(iPoint, bufDRecv[buf_offset+0], bufDRecv[buf_offset+1]); + if(nDim == 3) base_nodes->SetVelLapl_Z(iPoint, bufDRecv[buf_offset+2]); + break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index e149df25bd8..b65ee913686 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -190,6 +190,61 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain /*--- Upwind second order reconstruction and gradients ---*/ CommonPreprocessing(geometry, config, Output); + + if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED){ + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + + SU2_OMP_FOR_DYN(256) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + + const bool boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint); + + /*--- Initialize. ---*/ + nodes->SetVelLapl(iPoint, 0.0, 0.0); + if (nDim == 3) { + nodes->SetVelLapl_Z(iPoint, 0.0); + } + + /*--- Loop over the neighbors of point i. ---*/ + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + + bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint); + + /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ + if (boundary_i && !boundary_j) continue; + + /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ + const su2double distance = GeometryToolbox::Distance(nDim, nodes->GetMesh_Coord(iPoint), nodes->GetMesh_Coord(jPoint)); + + const su2double delta_x = (flowNodes->GetVelocity(jPoint,0)-flowNodes->GetVelocity(iPoint,0))/distance; + const su2double delta_y = (flowNodes->GetVelocity(jPoint,1)-flowNodes->GetVelocity(iPoint,1))/distance; + + nodes->AddVelLapl(iPoint, delta_x, delta_y); + if (nDim == 3) { + const su2double delta_z = (flowNodes->GetVelocity(jPoint,2)-flowNodes->GetVelocity(iPoint,2))/distance; + nodes->AddVelLapl_Z(iPoint, delta_z); + } + + } + + + } + END_SU2_OMP_FOR + + /*--- Correct the Laplacian across any periodic boundaries. ---*/ + + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VEL_LAPLACIAN); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VEL_LAPLACIAN); + } + + /*--- MPI parallelization ---*/ + + InitiateComms(geometry, config, VELOCITY_LAPLACIAN); + CompleteComms(geometry, config, VELOCITY_LAPLACIAN); + + } + } void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, @@ -364,12 +419,17 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); } + if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED){ + numerics->SetVelLapl(nodes->GetVelLapl_X(iPoint), nodes->GetVelLapl_Y(iPoint)); + if (nDim == 3) numerics->SetVelLapl_Z(nodes->GetVelLapl_Z(iPoint)); + } + /*--- Compute the source term ---*/ auto residual = numerics->ComputeResidual(config); /*--- Store the SAS function ---*/ - if (sstParsedOptions.sas) { + if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_SIMPLE) { nodes->SetFTrans(iPoint, numerics->GetFTrans()); } diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index 4a477665f23..a8ca6b18cf7 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -52,7 +52,12 @@ CTurbSSTVariable::CTurbSSTVariable(su2double kine, su2double omega, su2double mu muT.resize(nPoint) = mut; - if(sstParsedOptions.sas) FTrans.resize(nPoint) = su2double(1.0); + if(sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_SIMPLE) FTrans.resize(nPoint) = su2double(1.0); + if(sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED) { + VelocityLaplacian_X.resize(nPoint) = su2double(0.0); + VelocityLaplacian_Y.resize(nPoint) = su2double(0.0); + if (nDim == 3) VelocityLaplacian_Z.resize(nPoint) = su2double(0.0); + } } void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_viscosity, From a6be2cf6d44d88449e71a71b826c028f58458ef2 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 15 Sep 2023 17:53:45 +0200 Subject: [PATCH 03/16] - Fixed getting coordinates of points - Fixed DataType for Periodic Comms --- SU2_CFD/src/output/CFlowOutput.cpp | 4 ++-- SU2_CFD/src/solvers/CSolver.cpp | 11 ++++++++++- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 3 ++- 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index bf9d989d202..dc926357627 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1250,7 +1250,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy"); AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations"); - if (config->GetSSTParsedOptions().sas) AddVolumeOutput("FTRANS", "FTrans", "SOLUTION", "value of FTrans for SAS simulation"); + if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_SIMPLE) AddVolumeOutput("FTRANS", "FTrans", "SOLUTION", "value of FTrans for SAS simulation"); break; case TURB_FAMILY::NONE: @@ -1531,7 +1531,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, Node_Turb->GetSRSGridSize(iPoint)); - if (config->GetSSTParsedOptions().sas) SetVolumeOutputValue("FTRANS", iPoint, Node_Turb->GetFTrans(iPoint)); + if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_SIMPLE) SetVolumeOutputValue("FTRANS", iPoint, Node_Turb->GetFTrans(iPoint)); if (limiter) { SetVolumeOutputValue("LIMITER_TKE", iPoint, Node_Turb->GetLimiter(iPoint, 0)); SetVolumeOutputValue("LIMITER_DISSIPATION", iPoint, Node_Turb->GetLimiter(iPoint, 1)); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index a0afe86738e..dae8335567b 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -276,6 +276,11 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, MPI_TYPE = COMM_TYPE_DOUBLE; ICOUNT = nVar; break; + case PERIODIC_VEL_LAPLACIAN: + COUNT_PER_POINT = nDim; + MPI_TYPE = COMM_TYPE_DOUBLE; + ICOUNT = nDim; + break; default: SU2_MPI::Error("Unrecognized quantity for periodic communication.", CURRENT_FUNCTION); @@ -987,7 +992,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, if (!geometry->nodes->GetPeriodicBoundary(jPoint)) { /*--- Solution differences ---*/ - const su2double distance = GeometryToolbox::Distance(nDim, base_nodes->GetMesh_Coord(iPoint), base_nodes->GetMesh_Coord(jPoint)); + const su2double distance = GeometryToolbox::Distance(nDim, geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint)); for (iDim = 0; iDim < nDim; iDim++) Diff[iDim] = (base_nodes->GetVelocity(iPoint, iDim) - @@ -1442,6 +1447,10 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; break; + case VELOCITY_LAPLACIAN: + COUNT_PER_POINT = nDim; + MPI_TYPE = COMM_TYPE_DOUBLE; + break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index b65ee913686..da08a8f4903 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -28,6 +28,7 @@ #include "../../include/solvers/CTurbSSTSolver.hpp" #include "../../include/variables/CTurbSSTVariable.hpp" #include "../../include/variables/CFlowVariable.hpp" +#include "../../include/variables/CMeshVariable.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" @@ -214,7 +215,7 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain if (boundary_i && !boundary_j) continue; /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ - const su2double distance = GeometryToolbox::Distance(nDim, nodes->GetMesh_Coord(iPoint), nodes->GetMesh_Coord(jPoint)); + const su2double distance = GeometryToolbox::Distance(nDim, geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint)); const su2double delta_x = (flowNodes->GetVelocity(jPoint,0)-flowNodes->GetVelocity(iPoint,0))/distance; const su2double delta_y = (flowNodes->GetVelocity(jPoint,1)-flowNodes->GetVelocity(iPoint,1))/distance; From d1f4db63260f1eab2e71cda6c87352e160b0bdc1 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 19 Sep 2023 18:59:50 +0200 Subject: [PATCH 04/16] - Added output of laplacian --- SU2_CFD/src/output/CFlowOutput.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index dc926357627..2caa44a4ae6 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1251,6 +1251,11 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations"); if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_SIMPLE) AddVolumeOutput("FTRANS", "FTrans", "SOLUTION", "value of FTrans for SAS simulation"); + if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_COMPLICATED){ + AddVolumeOutput("VEL_LAPLACIAN_X", "Vel Laplacian x", "SOLUTION", "value of laplacian of velocity for SAS simulation"); + AddVolumeOutput("VEL_LAPLACIAN_Y", "Vel Laplacian y", "SOLUTION", "value of laplacian of velocity for SAS simulation"); + if (nDim == 3) AddVolumeOutput("VEL_LAPLACIAN_Z", "Vel Laplacian z", "SOLUTION", "value of laplacian of velocity for SAS simulation"); + } break; case TURB_FAMILY::NONE: @@ -1532,6 +1537,11 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, Node_Turb->GetSRSGridSize(iPoint)); if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_SIMPLE) SetVolumeOutputValue("FTRANS", iPoint, Node_Turb->GetFTrans(iPoint)); + if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_COMPLICATED){ + SetVolumeOutputValue("VEL_LAPLACIAN_X", iPoint, Node_Turb->GetVelLapl_X(iPoint)); + SetVolumeOutputValue("VEL_LAPLACIAN_Y", iPoint, Node_Turb->GetVelLapl_Y(iPoint)); + if (nDim == 3) SetVolumeOutputValue("VEL_LAPLACIAN_Z", iPoint, Node_Turb->GetVelLapl_Z(iPoint)); + } if (limiter) { SetVolumeOutputValue("LIMITER_TKE", iPoint, Node_Turb->GetLimiter(iPoint, 0)); SetVolumeOutputValue("LIMITER_DISSIPATION", iPoint, Node_Turb->GetLimiter(iPoint, 1)); From 35958b82fa9fa7345ea993c3250929bc89313649 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 19 Sep 2023 19:16:08 +0200 Subject: [PATCH 05/16] - Fixed missing variable initialization for CTurbSAVariable --- SU2_CFD/src/variables/CTurbSAVariable.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SU2_CFD/src/variables/CTurbSAVariable.cpp b/SU2_CFD/src/variables/CTurbSAVariable.cpp index 454ce5734a8..7e0e24756a0 100644 --- a/SU2_CFD/src/variables/CTurbSAVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSAVariable.cpp @@ -48,6 +48,9 @@ CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsi DES_LengthScale.resize(nPoint) = su2double(0.0); Vortex_Tilting.resize(nPoint); + + k.resize(nPoint) = su2double(0.0); + Omega.resize(nPoint) = su2double(0.0); } void CTurbSAVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView PrimGrad_Flow, From 09e977f35dd781bf6737d3b8e2a15c7bd14f440a Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 5 Oct 2023 18:13:05 +0200 Subject: [PATCH 06/16] - Added DES implementations to SST: DDES, Improved DDES (IDDES), Simplified IDDES (SIDDES) --- Common/include/option_structure.hpp | 8 +- SU2_CFD/include/numerics/CNumerics.hpp | 12 ++ .../numerics/turbulent/turb_sources.hpp | 3 + SU2_CFD/include/solvers/CTurbSSTSolver.hpp | 10 ++ SU2_CFD/include/variables/CTurbSAVariable.hpp | 14 --- .../include/variables/CTurbSSTVariable.hpp | 1 + SU2_CFD/include/variables/CTurbVariable.hpp | 14 +++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 118 ++++++++++++++++++ SU2_CFD/src/variables/CTurbSAVariable.cpp | 1 - SU2_CFD/src/variables/CTurbVariable.cpp | 2 + 10 files changed, 167 insertions(+), 16 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 5fe8228925e..96d7b55aa38 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1380,7 +1380,10 @@ enum ENUM_HYBRIDRANSLES { SA_DES = 1, /*!< \brief Kind of Hybrid RANS/LES (SA - Detached Eddy Simulation (DES)). */ SA_DDES = 2, /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Delta_max SGS ). */ SA_ZDES = 3, /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Vorticity based SGS like Zonal DES). */ - SA_EDDES = 4 /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES). */ + SA_EDDES = 4, /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES). */ + SST_DDES = 5, /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): DDES). */ + SST_IDDES = 6, /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Improved DDES). */ + SST_SIDDES = 7 /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Simplified Improved DDES). */ }; static const MapType HybridRANSLES_Map = { MakePair("NONE", NO_HYBRIDRANSLES) @@ -1388,6 +1391,9 @@ static const MapType HybridRANSLES_Map = { MakePair("SA_DDES", SA_DDES) MakePair("SA_ZDES", SA_ZDES) MakePair("SA_EDDES", SA_EDDES) + MakePair("SST_DDES", SST_DDES) + MakePair("SST_IDDES", SST_IDDES) + MakePair("SST_SIDDES", SST_SIDDES) }; /*! diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 14088f021cd..ffa7f1dec44 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 lengthScale_i, lengthScale_j; + public: /*! * \brief Return type used in some "ComputeResidual" overloads to give a @@ -828,6 +830,16 @@ class CNumerics { dist_j = val_dist_j; } + /*! + * \brief Set the value of the length scale for SST. + * \param[in] val_lengthScale_i - Value of of the length scale for SST from point i. + * \param[in] val_lengthScale_j - Value of of the length scale for SST from point j. + */ + void SetLengthScale(su2double val_lengthScale_i, su2double val_lengthScale_j) { + lengthScale_i = val_lengthScale_i; + lengthScale_j = val_lengthScale_j; + } + /*! * \brief Set the value of the roughness from the nearest wall. * \param[in] val_dist_i - Value of of the roughness of the nearest wall from point i diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index af920cd34c1..83d0ca9dfae 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -834,6 +834,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Dissipation ---*/ su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0]; + if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) + dk = Density_i * sqrt(ScalarVar_i[0]*ScalarVar_i[0]*ScalarVar_i[0]) / lengthScale_i; + su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1]; /*--- LM model coupling with production and dissipation term for k transport equation---*/ diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp index 6dbed05fd9d..aceaa94beef 100644 --- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp @@ -54,6 +54,16 @@ class CTurbSSTSolver final : public CTurbSolver { const CConfig *config, unsigned short val_marker); + /*! + * \brief A virtual member. + * \param[in] solver - Solver container + * \param[in] geometry - Geometrical definition. + * \param[in] config - Definition of the particular problem. + */ + void SetDES_LengthScale(CSolver** solver, + CGeometry *geometry, + CConfig *config); + public: /*! * \brief Constructor. diff --git a/SU2_CFD/include/variables/CTurbSAVariable.hpp b/SU2_CFD/include/variables/CTurbSAVariable.hpp index 569e416fa15..6a19003da02 100644 --- a/SU2_CFD/include/variables/CTurbSAVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp @@ -39,7 +39,6 @@ class CTurbSAVariable final : public CTurbVariable { private: - VectorType DES_LengthScale; VectorType Vortex_Tilting; public: @@ -60,19 +59,6 @@ class CTurbSAVariable final : public CTurbVariable { */ ~CTurbSAVariable() override = default; - /*! - * \brief Get the DES length scale - * \param[in] iPoint - Point index. - * \return Value of the DES length Scale. - */ - inline su2double GetDES_LengthScale(unsigned long iPoint) const override { return DES_LengthScale(iPoint); } - - /*! - * \brief Set the DES Length Scale. - * \param[in] iPoint - Point index. - */ - inline void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) override { DES_LengthScale(iPoint) = val_des_lengthscale; } - /*! * \brief Set the vortex tilting measure for computation of the EDDES length scale * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp index 0113a0ead3a..469679eb6bd 100644 --- a/SU2_CFD/include/variables/CTurbSSTVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSSTVariable.hpp @@ -87,4 +87,5 @@ class CTurbSSTVariable final : public CTurbVariable { * \brief Get the value of the cross diffusion of tke and omega. */ inline su2double GetCrossDiff(unsigned long iPoint) const override { return CDkw(iPoint); } + }; diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index 1a476fdac70..1196030b738 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -38,6 +38,7 @@ class CTurbVariable : public CScalarVariable { protected: VectorType muT; /*!< \brief Eddy viscosity. */ + VectorType DES_LengthScale; public: static constexpr size_t MAXNVAR = 2; @@ -100,5 +101,18 @@ class CTurbVariable : public CScalarVariable { */ inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; } + /*! + * \brief Get the DES length scale + * \param[in] iPoint - Point index. + * \return Value of the DES length Scale. + */ + inline su2double GetDES_LengthScale(unsigned long iPoint) const override { return DES_LengthScale(iPoint); } + + /*! + * \brief Set the DES Length Scale. + * \param[in] iPoint - Point index. + */ + inline void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) override { DES_LengthScale(iPoint) = val_des_lengthscale; } + }; diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index e1e9ef644c7..dc68a243763 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -188,8 +188,19 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);) + const auto kind_hybridRANSLES = config->GetKind_HybridRANSLES(); + /*--- Upwind second order reconstruction and gradients ---*/ CommonPreprocessing(geometry, config, Output); + + if (kind_hybridRANSLES != NO_HYBRIDRANSLES) { + + /*--- Compute the DES length scale ---*/ + + SetDES_LengthScale(solver_container, geometry, config); + + } + } void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, @@ -358,6 +369,10 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); } + if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) { + numerics->SetLengthScale(nodes->GetDES_LengthScale(iPoint), 0.0); + } + /*--- Compute the source term ---*/ auto residual = numerics->ComputeResidual(config); @@ -1060,3 +1075,106 @@ void CTurbSSTSolver::SetUniformInlet(const CConfig* config, unsigned short iMark } } + + +void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CConfig *config){ + + const auto kind_hybridRANSLES = config->GetKind_HybridRANSLES(); + + + auto* flowNodes = su2staticcast_p(solver[FLOW_SOL]->GetNodes()); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){ + const su2double StrainMag = max(flowNodes->GetStrainMag(iPoint), 1e-12); + const su2double VortMag = max(GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)), 1e-12); + + const su2double KolmConst2 = 0.41*0.41; + const su2double wallDist2 = geometry->nodes->GetWall_Distance(iPoint)*geometry->nodes->GetWall_Distance(iPoint); + + const su2double eddyVisc = nodes->GetmuT(iPoint)/flowNodes->GetDensity(iPoint); + const su2double lamVisc = nodes->GetLaminarViscosity(iPoint)/flowNodes->GetDensity(iPoint); + + const su2double C_DES1 = 0.78; + const su2double C_DES2 = 0.61; + + const su2double h_max = geometry->nodes->GetMaxLength(iPoint); + const su2double C_DES = C_DES1 * nodes->GetF1blending(iPoint) + C_DES2 * (1-nodes->GetF1blending(iPoint)); + const su2double l_RANS = sqrt(nodes->GetSolution(iPoint, 0)) / (constants[6] * nodes->GetSolution(iPoint, 1)); + + su2double DES_lengthScale = 0.0; + + switch(kind_hybridRANSLES){ + case SST_DDES: { + + const su2double r_d = (eddyVisc + lamVisc) / (KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))); + const su2double C_d1 = 20.0; + const su2double C_d2 = 3.0; + + const su2double f_d = 1 - tanh(pow(C_d1 * r_d, C_d2)); + + const su2double l_LES = C_DES * h_max; + + DES_lengthScale = l_RANS - f_d * max(0.0, l_RANS - l_LES); + + break; + } + case SST_IDDES: { + + // Constants + const su2double C_w = 0.15; + const su2double C_dt1 = 20.0; + const su2double C_dt2 = 3.0; + const su2double C_l = 5.0; + const su2double C_t = 1.87; + + const su2double alpha = 0.25 - sqrt(wallDist2) / h_max; + const su2double f_b = min(2.0 * exp(-9.0 * alpha*alpha), 1.0); + const su2double r_dt = eddyVisc / (KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))); + const su2double f_dt = 1 - tanh(pow(C_dt1 * r_dt, C_dt2)); + const su2double ftilda_d = max(1.0 - f_dt, f_b); + + const su2double r_dl = lamVisc / (KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))); + const su2double f_l = tanh(pow(C_l*C_l*r_dl, 10.0)); + const su2double f_t = tanh(pow(C_t*C_t*r_dt, 3.0)); + const su2double f_e2 = 1.0 - max(f_t, f_l); + const su2double f_e1 = alpha >= 0 ? 2.0 * exp(-11.09*alpha*alpha) : 2.0 * exp(-9.0*alpha*alpha); + const su2double f_e = f_e2 * max((f_e1 - 1.0), 0.0); + + + const su2double Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max); + const su2double l_LES = C_DES * Delta; + + DES_lengthScale = ftilda_d *(1.0+f_e)*l_RANS + (1.0 - ftilda_d) * l_LES; + + break; + } + case SST_SIDDES: { + + // Constants + const su2double C_w = 0.15; + const su2double C_dt1 = 20.0; + const su2double C_dt2 = 3.0; + const su2double C_l = 5.0; + const su2double C_t = 1.87; + + const su2double alpha = 0.25 - sqrt(wallDist2) / h_max; + const su2double f_b = min(2.0 * exp(-9.0 * alpha*alpha), 1.0); + const su2double r_dt = eddyVisc / (KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))); + const su2double f_dt = 1 - tanh(pow(C_dt1 * r_dt, C_dt2)); + const su2double ftilda_d = max(1.0 - f_dt, f_b); + + const su2double Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max); + const su2double l_LES = C_DES * Delta; + + DES_lengthScale = ftilda_d*l_RANS + (1.0 - ftilda_d) * l_LES; + + break; + } + } + + nodes->SetDES_LengthScale(iPoint, DES_lengthScale); + + } + END_SU2_OMP_FOR +} diff --git a/SU2_CFD/src/variables/CTurbSAVariable.cpp b/SU2_CFD/src/variables/CTurbSAVariable.cpp index 454ce5734a8..b8335969ca1 100644 --- a/SU2_CFD/src/variables/CTurbSAVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSAVariable.cpp @@ -46,7 +46,6 @@ CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsi Solution_time_n1 = Solution; } - DES_LengthScale.resize(nPoint) = su2double(0.0); Vortex_Tilting.resize(nPoint); } diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp index 54bede5033f..e905de5205a 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -35,4 +35,6 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned turb_index.resize(nPoint) = su2double(1.0); intermittency.resize(nPoint) = su2double(1.0); + DES_LengthScale.resize(nPoint) = su2double(0.0); + } From 8c8e89da6f5d6a36a48e5ab967bbdf4f4d226f17 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 5 Oct 2023 22:22:13 +0200 Subject: [PATCH 07/16] - Added SST DDES in config --- Common/src/CConfig.cpp | 3 +++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 2e8788fbe23..ffb186d0487 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6157,6 +6157,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case SA_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break; case SA_ZDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break; case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; + case SST_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; + case SST_IDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; + case SST_SIDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; } break; case MAIN_SOLVER::NEMO_EULER: diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index dc68a243763..871dde71555 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -1081,11 +1081,11 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C const auto kind_hybridRANSLES = config->GetKind_HybridRANSLES(); - auto* flowNodes = su2staticcast_p(solver[FLOW_SOL]->GetNodes()); SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){ + const su2double StrainMag = max(flowNodes->GetStrainMag(iPoint), 1e-12); const su2double VortMag = max(GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)), 1e-12); @@ -1138,7 +1138,7 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C const su2double f_l = tanh(pow(C_l*C_l*r_dl, 10.0)); const su2double f_t = tanh(pow(C_t*C_t*r_dt, 3.0)); const su2double f_e2 = 1.0 - max(f_t, f_l); - const su2double f_e1 = alpha >= 0 ? 2.0 * exp(-11.09*alpha*alpha) : 2.0 * exp(-9.0*alpha*alpha); + const su2double f_e1 = alpha >= 0.0 ? 2.0 * exp(-11.09*alpha*alpha) : 2.0 * exp(-9.0*alpha*alpha); const su2double f_e = f_e2 * max((f_e1 - 1.0), 0.0); From 71ac9c497a1117c6e2efdae2652ef3934a968a82 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 10 Oct 2023 18:01:55 +0200 Subject: [PATCH 08/16] - Added Debug quantities --- .../include/variables/CTurbSSTVariable.hpp | 39 +++++++++++++++++++ SU2_CFD/include/variables/CVariable.hpp | 19 +++++++++ SU2_CFD/src/output/CFlowOutput.cpp | 34 ++++++++++++++++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 6 +++ SU2_CFD/src/variables/CTurbSSTVariable.cpp | 7 ++++ 5 files changed, 105 insertions(+) diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp index 469679eb6bd..3e958314bb8 100644 --- a/SU2_CFD/include/variables/CTurbSSTVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSSTVariable.hpp @@ -45,6 +45,13 @@ class CTurbSSTVariable final : public CTurbVariable { VectorType F2; /*!< \brief Menter blending function for blending of k-w and k-eps. */ VectorType CDkw; /*!< \brief Cross-diffusion. */ SST_ParsedOptions sstParsedOptions; + + VectorType ftilda_d; + VectorType l_RANS; + VectorType l_LES; + VectorType r_dl; + VectorType r_dt; + VectorType r_d; public: /*! * \brief Constructor of the class. @@ -88,4 +95,36 @@ class CTurbSSTVariable final : public CTurbVariable { */ inline su2double GetCrossDiff(unsigned long iPoint) const override { return CDkw(iPoint); } + /*! + * \brief Set the DES Length Scale. + * \param[in] iPoint - Point index. + */ + inline void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_d) override { + ftilda_d(iPoint) = val_ftilda_d; + l_RANS(iPoint) = val_l_RANS; + l_LES(iPoint) = val_l_LES; + + if ( config->GetKind_HybridRANSLES() == SST_DDES) + r_d(iPoint) = val_r_d; + else + r_dt(iPoint) = val_r_d; + } + + inline void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_dl, su2double val_r_dt) override { + ftilda_d(iPoint) = val_ftilda_d; + l_RANS(iPoint) = val_l_RANS; + l_LES(iPoint) = val_l_LES; + r_dt(iPoint) = val_r_dt; + r_dl(iPoint) = val_r_dl; + } + + inline su2double Get_L_RANS(unsigned long iPoint) const override { return l_RANS(iPoint); } + inline su2double Get_L_LES(unsigned long iPoint) const override { return l_LES(iPoint); } + inline su2double Get_ftilda_d(unsigned long iPoint) const override { return ftilda_d(iPoint); } + inline su2double Get_r_dt(unsigned long iPoint) const override { return r_dt(iPoint); } + inline su2double Get_r_dl(unsigned long iPoint) const override { return r_dl(iPoint); } + inline su2double Get_r_d(unsigned long iPoint) const override { return r_d(iPoint); } + + + }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 050be2b7787..3c46655030a 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -381,6 +381,25 @@ class CVariable { */ inline virtual void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) {} + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + */ + inline virtual void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_d) {} + + inline virtual su2double Get_L_RANS(unsigned long iPoint) const { return 0.0; } + inline virtual su2double Get_L_LES(unsigned long iPoint) const { return 0.0; } + inline virtual su2double Get_ftilda_d(unsigned long iPoint) const { return 0.0; } + inline virtual su2double Get_r_dt(unsigned long iPoint) const { return 0.0; } + inline virtual su2double Get_r_dl(unsigned long iPoint) const { return 0.0; } + inline virtual su2double Get_r_d(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + */ + inline virtual void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_dl, su2double val_r_dt) {} + /*! * \brief A virtual member. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 75bd80feebb..5bfda2f6863 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1464,6 +1464,23 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) { if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) { AddVolumeOutput("DES_LENGTHSCALE", "DES_LengthScale", "DDES", "DES length scale value"); AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "DDES", "Wall distance value"); + if ( config->GetKind_HybridRANSLES() == SST_DDES){ + AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value"); + AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value"); + AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value"); + AddVolumeOutput("R_D", "r_d", "DDES", "DES length scale value"); + } else if ( config->GetKind_HybridRANSLES() == SST_IDDES){ + AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value"); + AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value"); + AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value"); + AddVolumeOutput("R_DT", "r_dt", "DDES", "DES length scale value"); + AddVolumeOutput("R_DL", "r_dl", "DDES", "DES length scale value"); + } else if ( config->GetKind_HybridRANSLES() == SST_SIDDES){ + AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value"); + AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value"); + AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value"); + AddVolumeOutput("R_DT", "r_dt", "DDES", "DES length scale value"); + } } if (config->GetViscous()) { @@ -1563,6 +1580,23 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) { SetVolumeOutputValue("DES_LENGTHSCALE", iPoint, Node_Flow->GetDES_LengthScale(iPoint)); SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint)); + if ( config->GetKind_HybridRANSLES() == SST_DDES){ + SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint)); + SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint)); + SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint)); + SetVolumeOutputValue("R_D", iPoint, Node_Turb->Get_r_d(iPoint)); + } else if ( config->GetKind_HybridRANSLES() == SST_IDDES){ + SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint)); + SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint)); + SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint)); + SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(iPoint)); + SetVolumeOutputValue("R_DL", iPoint, Node_Turb->Get_r_dl(iPoint)); + } else if ( config->GetKind_HybridRANSLES() == SST_SIDDES){ + SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint)); + SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint)); + SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint)); + SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(iPoint)); + } } switch (config->GetKind_Species_Model()) { diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 871dde71555..c99ebd53433 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -1117,6 +1117,8 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C DES_lengthScale = l_RANS - f_d * max(0.0, l_RANS - l_LES); + nodes->SetDebug_Quantities(config, iPoint, f_d, l_RANS, l_LES, r_d); + break; } case SST_IDDES: { @@ -1145,6 +1147,8 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C const su2double Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max); const su2double l_LES = C_DES * Delta; + nodes->SetDebug_Quantities(config, iPoint, ftilda_d, l_RANS, l_LES, r_dl, r_dt); + DES_lengthScale = ftilda_d *(1.0+f_e)*l_RANS + (1.0 - ftilda_d) * l_LES; break; @@ -1167,6 +1171,8 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C const su2double Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max); const su2double l_LES = C_DES * Delta; + nodes->SetDebug_Quantities(config, iPoint, ftilda_d, l_RANS, l_LES, r_dt); + DES_lengthScale = ftilda_d*l_RANS + (1.0 - ftilda_d) * l_LES; break; diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index 9778f1e37cc..959294d8643 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -50,6 +50,13 @@ CTurbSSTVariable::CTurbSSTVariable(su2double kine, su2double omega, su2double mu F2.resize(nPoint) = su2double(0.0); CDkw.resize(nPoint) = su2double(0.0); + ftilda_d.resize(nPoint) = su2double(0.0); + l_RANS.resize(nPoint) = su2double(0.0); + l_LES.resize(nPoint) = su2double(0.0); + r_dl.resize(nPoint) = su2double(0.0); + r_dt.resize(nPoint) = su2double(0.0); + r_d.resize(nPoint) = su2double(0.0); + muT.resize(nPoint) = mut; } From 5b5fc30ce150e006fcccef75bf94790f166ba682 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Tue, 17 Oct 2023 12:34:18 +0200 Subject: [PATCH 09/16] - Added division by 0 protection --- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index c99ebd53433..ca66c09cf3a 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -1107,7 +1107,7 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C switch(kind_hybridRANSLES){ case SST_DDES: { - const su2double r_d = (eddyVisc + lamVisc) / (KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))); + const su2double r_d = (eddyVisc + lamVisc) / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10); const su2double C_d1 = 20.0; const su2double C_d2 = 3.0; @@ -1132,11 +1132,11 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C const su2double alpha = 0.25 - sqrt(wallDist2) / h_max; const su2double f_b = min(2.0 * exp(-9.0 * alpha*alpha), 1.0); - const su2double r_dt = eddyVisc / (KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))); + const su2double r_dt = eddyVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10); const su2double f_dt = 1 - tanh(pow(C_dt1 * r_dt, C_dt2)); const su2double ftilda_d = max(1.0 - f_dt, f_b); - const su2double r_dl = lamVisc / (KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))); + const su2double r_dl = lamVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10); const su2double f_l = tanh(pow(C_l*C_l*r_dl, 10.0)); const su2double f_t = tanh(pow(C_t*C_t*r_dt, 3.0)); const su2double f_e2 = 1.0 - max(f_t, f_l); @@ -1164,7 +1164,7 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C const su2double alpha = 0.25 - sqrt(wallDist2) / h_max; const su2double f_b = min(2.0 * exp(-9.0 * alpha*alpha), 1.0); - const su2double r_dt = eddyVisc / (KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))); + const su2double r_dt = eddyVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10); const su2double f_dt = 1 - tanh(pow(C_dt1 * r_dt, C_dt2)); const su2double ftilda_d = max(1.0 - f_dt, f_b); From ef502830caebc7eb50892ad0212986dd742125a0 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sat, 21 Oct 2023 23:59:15 +0200 Subject: [PATCH 10/16] - Added computation of LESIQ to Output --- SU2_CFD/src/output/CFlowOutput.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 5bfda2f6863..780669ca2d8 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1481,6 +1481,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) { AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value"); AddVolumeOutput("R_DT", "r_dt", "DDES", "DES length scale value"); } + AddVolumeOutput("LESIQ", "LESIQ", "DDES", "LESIQ index for SRS simulations"); } if (config->GetViscous()) { @@ -1597,6 +1598,10 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint)); SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(iPoint)); } + const su2double mut = Node_Flow->GetEddyViscosity(iPoint); + const su2double mu = Node_Flow->GetLaminarViscosity(iPoint); + const su2double LESIQ = 1.0/(1.0+0.05*pow((mut+mu)/mu, 0.53)); + SetVolumeOutputValue("LESIQ", iPoint, LESIQ); } switch (config->GetKind_Species_Model()) { From 48a9f3dded8b00823ce6e6bd7b4e70e68bdcb83f Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 23 Oct 2023 19:07:58 +0200 Subject: [PATCH 11/16] - modified codi --- externals/codi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/externals/codi b/externals/codi index eee1b5eea2d..8ee822a9b0b 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit eee1b5eea2ded8126c34c1415e3b9cf15a3e70f2 +Subproject commit 8ee822a9b0bb8235a2494467b774e27fb64ff14f From ae33b055200ea5002e52a13912c043df6a25a797 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 27 Oct 2023 18:33:22 +0200 Subject: [PATCH 12/16] - Modified Velocity Laplacian (vector -> Matrix) - Changed names of SAS models - Added descriptions of functions - Clean up of output functions - Removed SA functions or SST variables --- Common/include/option_structure.hpp | 21 +++--- Common/src/CConfig.cpp | 17 ++++- SU2_CFD/include/numerics/CNumerics.hpp | 19 ++---- .../numerics/turbulent/turb_sources.hpp | 22 ++----- SU2_CFD/include/variables/CTurbSAVariable.hpp | 22 ------- .../include/variables/CTurbSSTVariable.hpp | 64 ++++++++++--------- SU2_CFD/include/variables/CVariable.hpp | 50 ++++----------- SU2_CFD/src/output/CFlowOutput.cpp | 42 ++++++------ SU2_CFD/src/solvers/CSolver.cpp | 19 +++--- SU2_CFD/src/solvers/CTurbSASolver.cpp | 26 -------- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 21 +++--- SU2_CFD/src/variables/CTurbSSTVariable.cpp | 8 +-- 12 files changed, 121 insertions(+), 210 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 4a459c947a4..4bd401df305 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -990,8 +990,8 @@ 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. */ - SAS_SIMPLE, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ - SAS_COMPLICATED, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ + SAS_TRAVIS, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ + SAS_BABU, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ }; static const MapType SST_Options_Map = { MakePair("NONE", SST_OPTIONS::NONE) @@ -1004,8 +1004,8 @@ static const MapType SST_Options_Map = { MakePair("VORTICITY", SST_OPTIONS::V) MakePair("KATO-LAUNDER", SST_OPTIONS::KL) MakePair("UQ", SST_OPTIONS::UQ) - MakePair("SAS_SIMPLE", SST_OPTIONS::SAS_SIMPLE) - MakePair("SAS_COMPLICATED", SST_OPTIONS::SAS_COMPLICATED) + MakePair("SAS_TRAVIS", SST_OPTIONS::SAS_TRAVIS) + MakePair("SAS_BABU", SST_OPTIONS::SAS_BABU) }; /*! @@ -1014,10 +1014,9 @@ static const MapType SST_Options_Map = { struct SST_ParsedOptions { SST_OPTIONS version = SST_OPTIONS::V1994; /*!< \brief Enum SST base model. */ SST_OPTIONS production = SST_OPTIONS::NONE; /*!< \brief Enum for production corrections/modifiers for SST model. */ - SST_OPTIONS sasModel = SST_OPTIONS::SAS_SIMPLE; /*!< \brief Enum SST base model. */ + SST_OPTIONS sasModel = SST_OPTIONS::NONE; /*!< \brief Enum SST base model. */ bool sust = false; /*!< \brief Bool for SST model with sustaining terms. */ bool uq = false; /*!< \brief Bool for using uncertainty quantification. */ - bool sas = false; /*!< \brief Bool for using Scale Adaptive Simulations. */ bool modified = false; /*!< \brief Bool for modified (m) SST model. */ }; @@ -1054,15 +1053,14 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne const bool sst_v = IsPresent(SST_OPTIONS::V); const bool sst_kl = IsPresent(SST_OPTIONS::KL); const bool sst_uq = IsPresent(SST_OPTIONS::UQ); - const bool sst_sas_simple = IsPresent(SST_OPTIONS::SAS_SIMPLE); - const bool sst_sas_comp = IsPresent(SST_OPTIONS::SAS_COMPLICATED); - const bool sst_sas = sst_sas_simple || sst_sas_comp; + const bool sst_sas_simple = IsPresent(SST_OPTIONS::SAS_TRAVIS); + const bool sst_sas_comp = IsPresent(SST_OPTIONS::SAS_BABU); if (sst_sas_simple && sst_sas_comp) { SU2_MPI::Error("Two versions (Simple and Complicated) selected for SAS under SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION); } else if (sst_sas_simple) { - SSTParsedOptions.sasModel = SST_OPTIONS::SAS_SIMPLE; + SSTParsedOptions.sasModel = SST_OPTIONS::SAS_TRAVIS; } else { - SSTParsedOptions.sasModel = SST_OPTIONS::SAS_COMPLICATED; + SSTParsedOptions.sasModel = SST_OPTIONS::SAS_BABU; } if (sst_1994 && sst_2003) { @@ -1087,7 +1085,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.sust = sst_sust; SSTParsedOptions.modified = sst_m; SSTParsedOptions.uq = sst_uq; - SSTParsedOptions.sas = sst_sas; return SSTParsedOptions; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index ffb186d0487..3b36d999754 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6116,6 +6116,17 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { break; } cout << "." << endl; + if(sstParsedOptions.sasModel != SST_OPTIONS::NONE) cout << "Scale Adaptive Simulation model: "; + switch (sstParsedOptions.sasModel) { + case SST_OPTIONS::SAS_TRAVIS: + cout << "Travis et al. (2004)"; + break; + case SST_OPTIONS::SAS_BABU: + cout << "Babu et al. (2016)"; + break; + default: + break; + } break; } switch (Kind_Trans_Model) { @@ -6157,9 +6168,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case SA_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break; case SA_ZDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break; case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; - case SST_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; - case SST_IDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; - case SST_SIDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; + case SST_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break; + case SST_IDDES: cout << "Improved Delayed Detached Eddy Simulation (IDDES) with Vorticity-based SGS" << endl; break; + case SST_SIDDES: cout << "Simplified Improved Delayed Detached Eddy Simulation (SIDDES) with Vorticity-based SGS" << endl; break; } break; case MAIN_SOLVER::NEMO_EULER: diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index e880ddfdaf0..995d8196706 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -136,6 +136,9 @@ class CNumerics { const su2double *TurbPsi_i, /*!< \brief Vector of adjoint turbulent variables at point i. */ *TurbPsi_j; /*!< \brief Vector of adjoint turbulent variables at point j. */ + su2double lengthScale_i, lengthScale_j; /*!< \brief length scale for SST */ + su2double FTrans; /*!< \brief SAS function */ + su2double *VelLapl; /*!< \brief Laplacian of the velocity */ CMatrixView ConsVar_Grad_i, /*!< \brief Gradient of conservative variables at point i. */ ConsVar_Grad_j, /*!< \brief Gradient of conservative variables at point j. */ @@ -190,10 +193,6 @@ class CNumerics { bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ - su2double lengthScale_i, lengthScale_j; - su2double FTrans; /*!< \brief SAS function */ - su2double VelLapl_X, VelLapl_Y, VelLapl_Z; - public: /*! * \brief Return type used in some "ComputeResidual" overloads to give a @@ -719,16 +718,8 @@ class CNumerics { /*! * \brief Get the value of the value of FTrans. */ - inline void SetVelLapl(su2double val_VelLapl_X, su2double val_VelLapl_Y) { - VelLapl_X = val_VelLapl_X; - VelLapl_Y = val_VelLapl_Y; - } - - /*! - * \brief Get the value of the value of FTrans. - */ - inline void SetVelLapl_Z(su2double val_VelLapl_Z) { - VelLapl_Z = val_VelLapl_Z; + inline void SetVelLapl(su2double* val_VelLapl) { + VelLapl = val_VelLapl; } /*! diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 584cb454574..3fcdf28be38 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -855,7 +855,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Compute SAS source term. ---*/ su2double Q_SAS = 0.0; - if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED) { + if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU) { const su2double KolmConst = 0.41; const su2double csi2 = 3.51; @@ -866,23 +866,13 @@ class CSourcePieceWise_TurbSST final : public CNumerics { // Scale of the modeled turbulence const su2double L = sqrt(ScalarVar_i[0]) / (pow(beta_star, 0.25) * ScalarVar_i[1]); // Von Karman Length Scale - su2double VelLaplMag = VelLapl_X*VelLapl_X + VelLapl_Y*VelLapl_Y; - if (nDim == 3) VelLaplMag += VelLapl_Z*VelLapl_Z; + const su2double VelLaplMag = GeometryToolbox::SquaredNorm(nDim, VelLapl); const su2double L_vK_1 = KolmConst * StrainMag_i / sqrt(VelLaplMag); const su2double L_vK_2 = C_s * sqrt(KolmConst * csi2 / (beta_blended/beta_star - alfa_blended)) * gridSize; const su2double L_vK = max(L_vK_1, L_vK_2); - su2double gradOmega = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) - gradOmega += ScalarVar_Grad_i[1][iDim]*ScalarVar_Grad_i[1][iDim]; - - gradOmega /= ScalarVar_i[1]; - - su2double gradTKE = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) - gradTKE += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim]; - - gradTKE /= ScalarVar_i[0]; + const su2double gradTKE = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[0])/ScalarVar_i[0]; + const su2double gradOmega = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[1])/ScalarVar_i[1]; const su2double Q_SAS_1 = csi2 * KolmConst * StrainMag_i * StrainMag_i * (L/L_vK) * (L/L_vK); const su2double Q_SAS_2 = C * (2*ScalarVar_i[0] / sigma_phi) * max(gradOmega, gradTKE); @@ -895,7 +885,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Add the dissipation terms to the residuals.---*/ FTrans = 1.0; - if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_SIMPLE) { + if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_TRAVIS) { FTrans = max(1.0, pow(StrainMag_i / (cTrans * VorticityMag), 2.0)); } Residual[0] -= dk * Volume * FTrans; @@ -916,7 +906,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume; - if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED) { + if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU) { Jacobian_i[0][0] += Q_SAS * Volume / ScalarVar_i[0]; } } diff --git a/SU2_CFD/include/variables/CTurbSAVariable.hpp b/SU2_CFD/include/variables/CTurbSAVariable.hpp index 953a9bb53a6..31c6e40f970 100644 --- a/SU2_CFD/include/variables/CTurbSAVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp @@ -75,26 +75,4 @@ class CTurbSAVariable final : public CTurbVariable { */ inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); } - /*! - * \brief Get the value of the turbulence kinetic energy. - * \return the value of the turbulence kinetic energy. - */ - inline su2double GetSSTVariables_k(unsigned long iPoint) const { return k(iPoint); } - - /*! - * \brief Get the value of the turbulence frequency Omega. - * \return the value of the turbulence frequency Omega. - */ - inline su2double GetSSTVariables_omega(unsigned long iPoint) const { return Omega(iPoint); } - - /*! - * \brief Set the value of the SST variables computed with SA solution. - * \param[in] val_k - * \param[in] val_Omega - */ - void SetSSTVariables(unsigned long iPoint, su2double val_k, su2double val_Omega) { - k(iPoint) = val_k; - Omega(iPoint) = val_Omega; - } - }; diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp index 5fba3f8d5f1..b2d9f0a1266 100644 --- a/SU2_CFD/include/variables/CTurbSSTVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSSTVariable.hpp @@ -53,9 +53,7 @@ class CTurbSSTVariable final : public CTurbVariable { VectorType r_dt; VectorType r_d; VectorType FTrans; - VectorType VelocityLaplacian_X; - VectorType VelocityLaplacian_Y; - VectorType VelocityLaplacian_Z; + MatrixType VelocityLaplacian; public: /*! @@ -133,59 +131,65 @@ class CTurbSSTVariable final : public CTurbVariable { /*! * \brief Get the value of the FTrans. + * \param[in] iPoint - Point index. */ inline su2double GetFTrans(unsigned long iPoint) const override { return FTrans(iPoint); } /*! * \brief Set the value of the FTrans. + * \param[in] iPoint - Point index. + * \param[in] val_FTrans - Value of the FTrans variable. */ inline void SetFTrans(unsigned long iPoint, su2double val_FTrans) override { FTrans(iPoint) = val_FTrans; }; /*! - * \brief Get the value of the FTrans. - */ - inline su2double GetVelLapl_X(unsigned long iPoint) const override { return VelocityLaplacian_X(iPoint); } - - /*! - * \brief Get the value of the FTrans. + * \brief Get the value of the velocity laplacian. + * \param[in] iPoint - Point index. */ - inline su2double GetVelLapl_Y(unsigned long iPoint) const override { return VelocityLaplacian_Y(iPoint); } + inline su2double* GetVelLapl(unsigned long iPoint) final { return VelocityLaplacian[iPoint]; } /*! - * \brief Get the value of the FTrans. + * \brief Get the value of the velocity laplacian. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. */ - inline su2double GetVelLapl_Z(unsigned long iPoint) const override { return VelocityLaplacian_Z(iPoint); } + inline su2double GetVelLapl(unsigned long iPoint, unsigned short iDim) const final { return VelocityLaplacian[iPoint][iDim]; } /*! - * \brief Set the value of the FTrans. - */ - inline void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) override { - VelocityLaplacian_X(iPoint) += val_VelLapl_X; - VelocityLaplacian_Y(iPoint) += val_VelLapl_Y; + * \brief Incrementally add the velocity laplacian vector. + * \param[in] iPoint - Point index. + * \param[in] val_VelLapl_X - X-Component of the velocity laplacian. + * \param[in] val_VelLapl_Y - Y-Component of the velocity laplacian. + * \param[in] val_VelLapl_Z - Z-Component of the velocity laplacian. + */ + inline void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y, su2double val_VelLapl_Z) override { + VelocityLaplacian(iPoint, 0) += val_VelLapl_X; + VelocityLaplacian(iPoint, 1) += val_VelLapl_Y; + if(nDim == 3) VelocityLaplacian(iPoint, 2) += val_VelLapl_Z; }; /*! - * \brief Set the value of the FTrans. + * \brief Set the value of the velocity laplacian. + * \param[in] iPoint - Point index. + * \param[in] val_VelLapl - Vector of velocity laplacian. */ - inline void AddVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) override { - VelocityLaplacian_Z(iPoint) += val_VelLapl_Z; + inline void SetVelLapl(unsigned long iPoint, su2double* val_VelLapl) override { + VelocityLaplacian(iPoint, 0) = val_VelLapl[0]; + VelocityLaplacian(iPoint, 1) = val_VelLapl[1]; + if(nDim == 3) VelocityLaplacian(iPoint, 2) = val_VelLapl[2]; }; /*! - * \brief Set the value of the FTrans. + * \brief Set the value of the velocity laplacian. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_VelLapl - value of the velocity laplacian. */ - inline void SetVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) override { - VelocityLaplacian_X(iPoint) = val_VelLapl_X; - VelocityLaplacian_Y(iPoint) = val_VelLapl_Y; + inline void SetVelLapl(unsigned long iPoint, unsigned short iDim, su2double val_VelLapl) override { + VelocityLaplacian(iPoint, iDim) = val_VelLapl; }; - /*! - * \brief Set the value of the FTrans. - */ - inline void SetVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) override { - VelocityLaplacian_Z(iPoint) = val_VelLapl_Z; - }; }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 6f784652d40..17e373f105e 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1677,34 +1677,29 @@ class CVariable { inline virtual void SetFTrans(unsigned long iPoint, su2double val_FTrans) {} /*! - * \brief Get the value of the value of FTrans. - */ - inline virtual su2double GetVelLapl_X(unsigned long iPoint) const { return 0.0; } - /*! - * \brief Get the value of the value of FTrans. - */ - inline virtual su2double GetVelLapl_Y(unsigned long iPoint) const { return 0.0; } - /*! - * \brief Get the value of the value of FTrans. + * \brief Get the value of the velocity laplacian. */ - inline virtual su2double GetVelLapl_Z(unsigned long iPoint) const { return 0.0; } + inline virtual su2double* GetVelLapl(unsigned long iPoint) { return nullptr; } /*! - * \brief Set the value of the value of FTrans. + * \brief Get the value of the velocity laplacian. */ - inline virtual void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) {} + inline virtual su2double GetVelLapl(unsigned long iPoint, unsigned short iDim) const { return 0.0; } + /*! - * \brief Set the value of the value of FTrans. + * \brief Incrementally add the velocity laplacian vector. */ - inline virtual void AddVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) {} + inline virtual void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y, su2double val_VelLapl_Z) {} + /*! - * \brief Set the value of the value of FTrans. + * \brief Set the value of the velocity laplacian. */ - inline virtual void SetVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) {} + inline virtual void SetVelLapl(unsigned long iPoint, su2double* val_VelLapl) {} + /*! - * \brief Set the value of the value of FTrans. + * \brief Set the value of the velocity laplacian. */ - inline virtual void SetVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) {} + inline virtual void SetVelLapl(unsigned long iPoint, unsigned short iDim, su2double val_VelLapl) {} /*! * \brief Get the value of the eddy viscosity. @@ -1754,25 +1749,6 @@ class CVariable { */ inline virtual void SetmuT(unsigned long iPoint, su2double val_muT) {} - /*! - * \brief Get the value of the turbulence kinetic energy. - * \return the value of the turbulence kinetic energy. - */ - inline virtual su2double GetSSTVariables_k(unsigned long iPoint) const { return 0.0; } - - /*! - * \brief Get the value of the Rate of dissipation Omega. - * \return the value of the Rate of dissipation Omega. - */ - inline virtual su2double GetSSTVariables_omega(unsigned long iPoint) const { return 0.0; } - - /*! - * \brief Set the value of the SST variables computed with SA solution. - * \param[in] val_k - * \param[in] val_Omega - */ - inline virtual void SetSSTVariables(unsigned long iPoint, su2double val_k, su2double val_Omega) {} - /*! * \brief Get the value of the desired grid size for Scale Resolving Simulations. * \return the value of the desired grid size for Scale Resolving Simulations. diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 5ca13276506..5b07264dbf7 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1250,11 +1250,11 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy"); AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations"); - if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_SIMPLE) AddVolumeOutput("FTRANS", "FTrans", "SOLUTION", "value of FTrans for SAS simulation"); - if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_COMPLICATED){ - AddVolumeOutput("VEL_LAPLACIAN_X", "Vel Laplacian x", "SOLUTION", "value of laplacian of velocity for SAS simulation"); - AddVolumeOutput("VEL_LAPLACIAN_Y", "Vel Laplacian y", "SOLUTION", "value of laplacian of velocity for SAS simulation"); - if (nDim == 3) AddVolumeOutput("VEL_LAPLACIAN_Z", "Vel Laplacian z", "SOLUTION", "value of laplacian of velocity for SAS simulation"); + if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_TRAVIS) AddVolumeOutput("FTRANS", "FTrans", "SOLUTION", "value of FTrans for SAS simulation"); + if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_BABU){ + AddVolumeOutput("VEL-LAPLACIAN_X", "Vel Laplacian x", "SOLUTION", "value of laplacian of x-velocity for SAS simulation"); + AddVolumeOutput("VEL-LAPLACIAN_Y", "Vel Laplacian y", "SOLUTION", "value of laplacian of y-velocity for SAS simulation"); + if (nDim == 3) AddVolumeOutput("VEL-LAPLACIAN_Z", "Vel Laplacian z", "SOLUTION", "value of laplacian of z-velocity for SAS simulation"); } break; @@ -1472,22 +1472,18 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) { if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) { AddVolumeOutput("DES_LENGTHSCALE", "DES_LengthScale", "DDES", "DES length scale value"); AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "DDES", "Wall distance value"); + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + AddVolumeOutput("F_D", "f_d", "DDES", "Empiric blending function"); + AddVolumeOutput("L_RANS", "l_RANS", "DDES", "RANS length scale value"); + AddVolumeOutput("L_LES", "l_LES", "DDES", "LES length scale value"); + } if ( config->GetKind_HybridRANSLES() == SST_DDES){ - AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value"); - AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value"); - AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value"); - AddVolumeOutput("R_D", "r_d", "DDES", "DES length scale value"); + AddVolumeOutput("R_D", "r_d", "DDES", "r_d"); } else if ( config->GetKind_HybridRANSLES() == SST_IDDES){ - AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value"); - AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value"); - AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value"); - AddVolumeOutput("R_DT", "r_dt", "DDES", "DES length scale value"); - AddVolumeOutput("R_DL", "r_dl", "DDES", "DES length scale value"); + AddVolumeOutput("R_DT", "r_dt", "DDES", "turbulent r_d"); + AddVolumeOutput("R_DL", "r_dl", "DDES", "laminar r_d"); } else if ( config->GetKind_HybridRANSLES() == SST_SIDDES){ - AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value"); - AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value"); - AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value"); - AddVolumeOutput("R_DT", "r_dt", "DDES", "DES length scale value"); + AddVolumeOutput("R_DT", "r_dt", "DDES", "turbulent r_d"); } AddVolumeOutput("LESIQ", "LESIQ", "DDES", "LESIQ index for SRS simulations"); } @@ -1554,11 +1550,11 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, Node_Turb->GetSRSGridSize(iPoint)); - if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_SIMPLE) SetVolumeOutputValue("FTRANS", iPoint, Node_Turb->GetFTrans(iPoint)); - if (config->GetSSTParsedOptions().sas && config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_COMPLICATED){ - SetVolumeOutputValue("VEL_LAPLACIAN_X", iPoint, Node_Turb->GetVelLapl_X(iPoint)); - SetVolumeOutputValue("VEL_LAPLACIAN_Y", iPoint, Node_Turb->GetVelLapl_Y(iPoint)); - if (nDim == 3) SetVolumeOutputValue("VEL_LAPLACIAN_Z", iPoint, Node_Turb->GetVelLapl_Z(iPoint)); + if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_TRAVIS) SetVolumeOutputValue("FTRANS", iPoint, Node_Turb->GetFTrans(iPoint)); + if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_BABU){ + SetVolumeOutputValue("VEL-LAPLACIAN_X", iPoint, Node_Turb->GetVelLapl(iPoint, 0)); + SetVolumeOutputValue("VEL-LAPLACIAN_Y", iPoint, Node_Turb->GetVelLapl(iPoint, 1)); + if (nDim == 3) SetVolumeOutputValue("VEL-LAPLACIAN_Z", iPoint, Node_Turb->GetVelLapl(iPoint, 2)); } if (limiter) { SetVolumeOutputValue("LIMITER_TKE", iPoint, Node_Turb->GetLimiter(iPoint, 0)); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index dae8335567b..b1c10838a40 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1350,15 +1350,15 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, break; case PERIODIC_VEL_LAPLACIAN: - + { /*--- Adjust the undivided Laplacian. The accumulation was with a subtraction before communicating, so now just add. ---*/ - - base_nodes->AddVelLapl(iPoint, bufDRecv[buf_offset+0], bufDRecv[buf_offset+1]); - if(nDim == 3) base_nodes->AddVelLapl_Z(iPoint, bufDRecv[buf_offset+2]); + su2double bufDRecv_Z = 0.0; + if(nDim == 3) bufDRecv_Z = bufDRecv[buf_offset+2]; + base_nodes->AddVelLapl(iPoint, bufDRecv[buf_offset+0], bufDRecv[buf_offset+1], bufDRecv_Z); break; - + } default: SU2_MPI::Error("Unrecognized quantity for periodic communication.", @@ -1601,9 +1601,8 @@ void CSolver::InitiateComms(CGeometry *geometry, bufDSend[buf_offset+iVar] = base_nodes->GetSolution_time_n1(iPoint, iVar); break; case VELOCITY_LAPLACIAN: - bufDSend[buf_offset+0] = base_nodes->GetVelLapl_X(iPoint); - bufDSend[buf_offset+1] = base_nodes->GetVelLapl_Y(iPoint); - if(nDim == 3) bufDSend[buf_offset+2] = base_nodes->GetVelLapl_Z(iPoint); + for(iDim = 0; iDim < nDim; iDim++) + bufDSend[buf_offset+iDim] = base_nodes->GetVelLapl(iPoint, iDim); break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", @@ -1754,8 +1753,8 @@ void CSolver::CompleteComms(CGeometry *geometry, base_nodes->Set_Solution_time_n1(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; case VELOCITY_LAPLACIAN: - base_nodes->SetVelLapl(iPoint, bufDRecv[buf_offset+0], bufDRecv[buf_offset+1]); - if(nDim == 3) base_nodes->SetVelLapl_Z(iPoint, bufDRecv[buf_offset+2]); + for (iDim = 0; iDim < nDim; iDim++) + base_nodes->SetVelLapl(iPoint, iDim, bufDRecv[buf_offset+iDim]); break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index a81f75a5e09..65b56ab5d3d 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -240,32 +240,6 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain nodes->SetmuT(iPoint,muT); - // Now compute k and w - auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); - - const su2double BetaStar = 0.09; - const su2double a1 = 0.31; - - const su2double Omega = GeometryToolbox::Norm(nDim, flowNodes->GetVorticity(iPoint))/sqrt(BetaStar); - - const su2double nut = muT / rho; - const su2double arg2_1 = 2.0 * sqrt(nut / Omega) / (BetaStar * wallDistance); - const su2double arg2_2 = 500.0 * flowNodes->GetLaminarViscosity(iPoint) / (wallDistance*wallDistance * Omega); - - const su2double arg2 = max(arg2_1, arg2_2); - const su2double F2 = tanh(arg2*arg2); - - const su2double k = nut * max(Omega, flowNodes->GetStrainMag(iPoint) * F2 / a1); - - nodes->SetSSTVariables(iPoint, k, Omega); - - // Now compute desired cell size for Scale Resolving Simulations - const su2double Cmu = 0.09; - const su2double RANSLength = sqrt(k) / max(1e-20, (Cmu * Omega)); - const su2double RatioL = 0.1; // it should be less or equal than 0.2 - 0.1. Should be taken as input from config? - const su2double SRSGridSize = RANSLength * RatioL; - nodes->SetSRSGridSize(iPoint, SRSGridSize); - } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index ca33f730f22..c1a146f65f6 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -202,7 +202,7 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain } - if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED){ + if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU){ auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); SU2_OMP_FOR_DYN(256) @@ -211,10 +211,8 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain const bool boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint); /*--- Initialize. ---*/ - nodes->SetVelLapl(iPoint, 0.0, 0.0); - if (nDim == 3) { - nodes->SetVelLapl_Z(iPoint, 0.0); - } + for (unsigned short iDim = 0; iDim < nDim; iDim++) + nodes->SetVelLapl(iPoint, iDim, 0.0); /*--- Loop over the neighbors of point i. ---*/ for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { @@ -229,12 +227,12 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain const su2double delta_x = (flowNodes->GetVelocity(jPoint,0)-flowNodes->GetVelocity(iPoint,0))/distance; const su2double delta_y = (flowNodes->GetVelocity(jPoint,1)-flowNodes->GetVelocity(iPoint,1))/distance; + su2double delta_z = 0.0; - nodes->AddVelLapl(iPoint, delta_x, delta_y); if (nDim == 3) { - const su2double delta_z = (flowNodes->GetVelocity(jPoint,2)-flowNodes->GetVelocity(iPoint,2))/distance; - nodes->AddVelLapl_Z(iPoint, delta_z); + delta_z = (flowNodes->GetVelocity(jPoint,2)-flowNodes->GetVelocity(iPoint,2))/distance; } + nodes->AddVelLapl(iPoint, delta_x, delta_y, delta_z); } @@ -434,9 +432,8 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta numerics->SetLengthScale(nodes->GetDES_LengthScale(iPoint), 0.0); } - if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED){ - numerics->SetVelLapl(nodes->GetVelLapl_X(iPoint), nodes->GetVelLapl_Y(iPoint)); - if (nDim == 3) numerics->SetVelLapl_Z(nodes->GetVelLapl_Z(iPoint)); + if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU){ + numerics->SetVelLapl(nodes->GetVelLapl(iPoint)); } /*--- Compute the source term ---*/ @@ -444,7 +441,7 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta auto residual = numerics->ComputeResidual(config); /*--- Store the SAS function ---*/ - if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_SIMPLE) { + if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_TRAVIS) { nodes->SetFTrans(iPoint, numerics->GetFTrans()); } diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index f3455e4f1aa..e3b38e14e61 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -59,11 +59,9 @@ CTurbSSTVariable::CTurbSSTVariable(su2double kine, su2double omega, su2double mu muT.resize(nPoint) = mut; - if(sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_SIMPLE) FTrans.resize(nPoint) = su2double(1.0); - if(sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED) { - VelocityLaplacian_X.resize(nPoint) = su2double(0.0); - VelocityLaplacian_Y.resize(nPoint) = su2double(0.0); - if (nDim == 3) VelocityLaplacian_Z.resize(nPoint) = su2double(0.0); + if(sstParsedOptions.sasModel == SST_OPTIONS::SAS_TRAVIS) FTrans.resize(nPoint) = su2double(1.0); + if(sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU) { + VelocityLaplacian.resize(nPoint, nDim) = su2double(0.0); } } From 6443070d2f31f6e3a976204f2c54c6edb4971c98 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 20 Sep 2024 09:37:59 +0200 Subject: [PATCH 13/16] - Added outputs log for SST-DDES --- Common/src/CConfig.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 3b36d999754..e7f4500b41d 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6168,9 +6168,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case SA_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break; case SA_ZDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break; case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; - case SST_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break; - case SST_IDDES: cout << "Improved Delayed Detached Eddy Simulation (IDDES) with Vorticity-based SGS" << endl; break; - case SST_SIDDES: cout << "Simplified Improved Delayed Detached Eddy Simulation (SIDDES) with Vorticity-based SGS" << endl; break; + case SST_DDES: cout << "Delayed Detached Eddy Simulation (DDES)" << endl; break; + case SST_IDDES: cout << "Improved Delayed Detached Eddy Simulation (IDDES)" << endl; break; + case SST_SIDDES: cout << "Simplified Improved Delayed Detached Eddy Simulation (SIDDES)" << endl; break; } break; case MAIN_SOLVER::NEMO_EULER: From a58452681184e589853b4a0086de5a6cd948e991 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Fri, 20 Sep 2024 14:46:57 +0200 Subject: [PATCH 14/16] - update externals --- externals/medi | 2 +- externals/opdi | 2 +- subprojects/MLPCpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/externals/medi b/externals/medi index ab3a7688f6d..aafc2d1966b 160000 --- a/externals/medi +++ b/externals/medi @@ -1 +1 @@ -Subproject commit ab3a7688f6d518f8d940eb61a341d89f51922ba4 +Subproject commit aafc2d1966ba1233640af737e71c77c1a86183fd diff --git a/externals/opdi b/externals/opdi index 8c897988172..c42cca71a3d 160000 --- a/externals/opdi +++ b/externals/opdi @@ -1 +1 @@ -Subproject commit 8c89798817253abb017d857a0ae7f0520187645c +Subproject commit c42cca71a3d0b44fb482e268ecd40b623e71776b diff --git a/subprojects/MLPCpp b/subprojects/MLPCpp index c19c53ea2b8..665c45b7d35 160000 --- a/subprojects/MLPCpp +++ b/subprojects/MLPCpp @@ -1 +1 @@ -Subproject commit c19c53ea2b85ccfb185f1c6c87044dc0b5bc7ae0 +Subproject commit 665c45b7d3533c977eb1f637918d5b8b75c07d3b From 710e1dc5c1ec58128368db805b7348c3d12fdee2 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Oct 2024 14:59:26 +0100 Subject: [PATCH 15/16] Fix SAS Babu implementation --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 0204ad3b774..8b2bca6bcbd 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -940,8 +940,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics { const su2double L_vK_2 = C_s * sqrt(KolmConst * csi2 / (beta_blended/beta_star - alfa_blended)) * gridSize; const su2double L_vK = max(L_vK_1, L_vK_2); - const su2double gradTKE = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[0])/ScalarVar_i[0]; - const su2double gradOmega = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[1])/ScalarVar_i[1]; + const su2double gradTKE = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[0])/(ScalarVar_i[0]*ScalarVar_i[0]); + const su2double gradOmega = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[1])/(ScalarVar_i[1]*ScalarVar_i[1]); const su2double Q_SAS_1 = csi2 * KolmConst * StrainMag_i * StrainMag_i * (L/L_vK) * (L/L_vK); const su2double Q_SAS_2 = C * (2*ScalarVar_i[0] / sigma_phi) * max(gradOmega, gradTKE); From e4f28d4615f39f300da4e4672bedadfd545c9146 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Oct 2024 15:38:20 +0100 Subject: [PATCH 16/16] - Fixed compilation error --- SU2_CFD/src/solvers/CSolver.cpp | 6 +++--- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index c3550371a45..d19df52ab70 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1446,7 +1446,7 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case VELOCITY_LAPLACIAN: + case MPI_QUANTITIES::VELOCITY_LAPLACIAN: COUNT_PER_POINT = nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; @@ -1599,7 +1599,7 @@ void CSolver::InitiateComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution_time_n1(iPoint, iVar); break; - case VELOCITY_LAPLACIAN: + case MPI_QUANTITIES::VELOCITY_LAPLACIAN: for(iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset+iDim] = base_nodes->GetVelLapl(iPoint, iDim); break; @@ -1751,7 +1751,7 @@ void CSolver::CompleteComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) base_nodes->Set_Solution_time_n1(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; - case VELOCITY_LAPLACIAN: + case MPI_QUANTITIES::VELOCITY_LAPLACIAN: for (iDim = 0; iDim < nDim; iDim++) base_nodes->SetVelLapl(iPoint, iDim, bufDRecv[buf_offset+iDim]); break; diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index ac6c82a184c..20f5f5602e2 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -258,8 +258,8 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain /*--- MPI parallelization ---*/ - InitiateComms(geometry, config, VELOCITY_LAPLACIAN); - CompleteComms(geometry, config, VELOCITY_LAPLACIAN); + InitiateComms(geometry, config, MPI_QUANTITIES::VELOCITY_LAPLACIAN); + CompleteComms(geometry, config, MPI_QUANTITIES::VELOCITY_LAPLACIAN); }