Skip to content

Commit 0c8d49a

Browse files
author
rois1995
committed
- Fixed some notation
- Used GeometryToolbox when possible
1 parent ac2078c commit 0c8d49a

File tree

8 files changed

+68
-66
lines changed

8 files changed

+68
-66
lines changed

‎Common/include/option_structure.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1480,15 +1480,15 @@ enum ENUM_HYBRIDRANSLES {
14801480
SST_IDDES = 6, /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Improved DDES). */
14811481
SST_SIDDES = 7, /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Simplified Improved DDES). */
14821482
SST_EDDES = 8, /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
1483-
SA_EDDES_MOD = 9 /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
1483+
SA_EDDES_UNSTR = 9 /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
14841484
};
14851485
static const MapType<std::string, ENUM_HYBRIDRANSLES> HybridRANSLES_Map = {
14861486
MakePair("NONE", NO_HYBRIDRANSLES)
14871487
MakePair("SA_DES", SA_DES)
14881488
MakePair("SA_DDES", SA_DDES)
14891489
MakePair("SA_ZDES", SA_ZDES)
14901490
MakePair("SA_EDDES", SA_EDDES)
1491-
MakePair("SA_EDDES_MOD", SA_EDDES_MOD)
1491+
MakePair("SA_EDDES_UNSTR", SA_EDDES_UNSTR)
14921492
MakePair("SST_DDES", SST_DDES)
14931493
MakePair("SST_IDDES", SST_IDDES)
14941494
MakePair("SST_SIDDES", SST_SIDDES)

‎Common/src/CConfig.cpp‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -6442,13 +6442,14 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
64426442
cout << "Hybrid RANS/LES: ";
64436443
switch (Kind_HybridRANSLES) {
64446444
case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break;
6445-
case SA_DES: cout << "Detached Eddy Simulation (DES97) " << endl; break;
6446-
case SA_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
6447-
case SA_ZDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
6448-
case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
6449-
case SST_DDES: cout << "Delayed Detached Eddy Simulation (DDES)" << endl; break;
6450-
case SST_IDDES: cout << "Improved Delayed Detached Eddy Simulation (IDDES)" << endl; break;
6451-
case SST_SIDDES: cout << "Simplified Improved Delayed Detached Eddy Simulation (SIDDES)" << endl; break;
6445+
case SA_DES: cout << "Detached Eddy Simulation (DES97) " << endl; break;
6446+
case SA_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
6447+
case SA_ZDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
6448+
case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
6449+
case SA_EDDES_UNSTR: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS (modified for Unstructured grids)" << endl; break;
6450+
case SST_DDES: cout << "Delayed Detached Eddy Simulation (DDES)" << endl; break;
6451+
case SST_IDDES: cout << "Improved Delayed Detached Eddy Simulation (IDDES)" << endl; break;
6452+
case SST_SIDDES: cout << "Simplified Improved Delayed Detached Eddy Simulation (SIDDES)" << endl; break;
64526453
}
64536454
break;
64546455
case MAIN_SOLVER::NEMO_EULER:

‎SU2_CFD/include/numerics/turbulent/turb_sources.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -887,17 +887,17 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
887887
}
888888

889889
if (sstParsedOptions.production == SST_OPTIONS::COMP_Sarkar) {
890-
const su2double Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * Mt * Mt;
890+
const su2double Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * pow(Mt, 2);
891891
pk += Dilatation_Sarkar;
892892
}
893893

894894
/*--- Dissipation ---*/
895895

896896
su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt);
897897
if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES)
898-
dk = Density_i * sqrt(ScalarVar_i[0]*ScalarVar_i[0]*ScalarVar_i[0]) / lengthScale_i;
898+
dk = Density_i * sqrt(pow(ScalarVar_i[0], 3)) / lengthScale_i;
899899

900-
su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - 0.09/beta_blended * zetaFMt);
900+
su2double dw = beta_blended * Density_i * pow(ScalarVar_i[1], 2) * (1.0 - 0.09/beta_blended * zetaFMt);
901901

902902
/*--- LM model coupling with production and dissipation term for k transport equation---*/
903903
if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) {

‎SU2_CFD/include/variables/CTurbVariable.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -134,7 +134,7 @@ class CTurbVariable : public CScalarVariable {
134134
* \brief Set the vortex tilting measure for computation of the EDDES length scale
135135
* \param[in] iPoint - Point index.
136136
*/
137-
void SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
137+
void SetVortex_Tilting(unsigned long iPoint, const su2double Strain[3][3],
138138
const su2double* Vorticity, su2double LaminarViscosity) override;
139139

140140
/*!

‎SU2_CFD/include/variables/CVariable.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2247,7 +2247,7 @@ class CVariable {
22472247

22482248
inline virtual su2double GetTau_Wall(unsigned long iPoint) const { return 0.0; }
22492249

2250-
inline virtual void SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
2250+
inline virtual void SetVortex_Tilting(unsigned long iPoint, const su2double Strain[3][3],
22512251
const su2double* Vorticity, su2double LaminarViscosity) {}
22522252

22532253
inline virtual su2double GetVortex_Tilting(unsigned long iPoint) const { return 0.0; }

‎SU2_CFD/src/solvers/CTurbSASolver.cpp‎

Lines changed: 22 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -186,15 +186,21 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe
186186

187187
/*--- Set the vortex tilting coefficient at every node if required ---*/
188188

189-
if (kind_hybridRANSLES == SA_EDDES || kind_hybridRANSLES == SA_EDDES_MOD){
189+
if (kind_hybridRANSLES == SA_EDDES || kind_hybridRANSLES == SA_EDDES_UNSTR){
190190
auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
191191

192192
SU2_OMP_FOR_STAT(omp_chunk_size)
193193
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
194+
su2double Grad_Vel[3][3] = {{0.0}}, StrainMat[3][3] = {{0.0}};
194195
auto Vorticity = flowNodes->GetVorticity(iPoint);
195-
auto PrimGrad_Flow = flowNodes->GetGradient_Primitive(iPoint);
196+
for (unsigned short iDim = 0; iDim < nDim; iDim++) {
197+
for (unsigned short jDim = 0; jDim < nDim; jDim++) {
198+
Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Velocity() + iDim, jDim);
199+
}
200+
}
196201
auto Laminar_Viscosity = flowNodes->GetLaminarViscosity(iPoint);
197-
nodes->SetVortex_Tilting(iPoint, PrimGrad_Flow, Vorticity, Laminar_Viscosity);
202+
CNumerics::ComputeMeanRateOfStrainMatrix(3, StrainMat, Grad_Vel);
203+
nodes->SetVortex_Tilting(iPoint, StrainMat, Vorticity, Laminar_Viscosity);
198204
}
199205
END_SU2_OMP_FOR
200206
}
@@ -1456,8 +1462,8 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
14561462
pow(ratioOmega[1], 2)*delta[0]*delta[2] +
14571463
pow(ratioOmega[2], 2)*delta[0]*delta[1]);
14581464

1459-
const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0));
1460-
const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0));
1465+
const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2));
1466+
const su2double f_d = 1.0-tanh(pow(8.0*r_d,3));
14611467

14621468
if (f_d < 0.99){
14631469
maxDelta = deltaDDES;
@@ -1507,8 +1513,8 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
15071513
min(f_max,
15081514
f_min + ((f_max - f_min)/(a2 - a1)) * (vortexTiltingMeasure - a1)));
15091515

1510-
const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0));
1511-
const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0));
1516+
const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2));
1517+
const su2double f_d = 1.0-tanh(pow(8.0*r_d,3));
15121518

15131519
su2double maxDelta = (ln_max/sqrt(3.0)) * f_kh;
15141520
if (f_d < 0.999){
@@ -1520,7 +1526,7 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
15201526

15211527
break;
15221528
}
1523-
case SA_EDDES_MOD: {
1529+
case SA_EDDES_UNSTR: {
15241530
/*--- An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows.
15251531
Shur et al.
15261532
Flow Turbulence Combust - 2015
@@ -1531,13 +1537,13 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
15311537
const su2double omega = GeometryToolbox::Norm(3, vorticity);
15321538

15331539
su2double ratioOmega[MAXNDIM] = {};
1534-
1535-
for (auto iDim = 0; iDim < 3; iDim++){
1540+
for (auto iDim = 0; iDim < MAXNDIM; iDim++){
15361541
ratioOmega[iDim] = vorticity[iDim]/omega;
15371542
}
15381543

15391544
const su2double deltaDDES = geometry->nodes->GetMaxLength(iPoint);
1540-
1545+
1546+
// Introduced correction for unstructured grids
15411547
su2double ln_max = -1;
15421548
for (const auto jPoint : geometry->nodes->GetPoints(iPoint)){
15431549
const auto coord_j = geometry->nodes->GetCoord(jPoint);
@@ -1546,13 +1552,13 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
15461552
const auto coord_k = geometry->nodes->GetCoord(kPoint);
15471553

15481554
su2double delta[MAXNDIM] = {};
1549-
// This should only be performed on 3D cases anyway
1550-
for (auto iDim = 0u; iDim < 3; iDim++){
1551-
delta[iDim] = (coord_j[iDim] - coord_k[iDim])/2.0; // Should I divide by 2 as I am interested in the dual volume?
1555+
for (auto iDim = 0u; iDim < MAXNDIM; iDim++){
1556+
// TODO: Should I divide by 2 as I am interested in the dual volume (the edge is split at midpoint)?
1557+
delta[iDim] = (coord_j[iDim] - coord_k[iDim])/2.0;
15521558
}
1553-
su2double l_n_minus_m[3];
1559+
su2double l_n_minus_m[MAXNDIM];
15541560
GeometryToolbox::CrossProduct(delta, ratioOmega, l_n_minus_m);
1555-
ln_max = max(ln_max, GeometryToolbox::Norm(nDim, l_n_minus_m));
1561+
ln_max = max(ln_max, GeometryToolbox::Norm(3, l_n_minus_m));
15561562
}
15571563
vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
15581564
}

‎SU2_CFD/src/solvers/CTurbSSTSolver.cpp‎

Lines changed: 17 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,6 @@
2828
#include "../../include/solvers/CTurbSSTSolver.hpp"
2929
#include "../../include/variables/CTurbSSTVariable.hpp"
3030
#include "../../include/variables/CFlowVariable.hpp"
31-
#include "../../include/variables/CMeshVariable.hpp"
3231
#include "../../../Common/include/parallelization/omp_structure.hpp"
3332
#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
3433

@@ -214,10 +213,16 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain
214213

215214
SU2_OMP_FOR_STAT(omp_chunk_size)
216215
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
216+
su2double Grad_Vel[3][3] = {{0.0}}, StrainMat[3][3] = {{0.0}};
217217
auto Vorticity = flowNodes->GetVorticity(iPoint);
218-
auto PrimGrad_Flow = flowNodes->GetGradient_Primitive(iPoint);
218+
for (unsigned short iDim = 0; iDim < nDim; iDim++) {
219+
for (unsigned short jDim = 0; jDim < nDim; jDim++) {
220+
Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Velocity() + iDim, jDim);
221+
}
222+
}
219223
auto Laminar_Viscosity = flowNodes->GetLaminarViscosity(iPoint);
220-
nodes->SetVortex_Tilting(iPoint, PrimGrad_Flow, Vorticity, Laminar_Viscosity);
224+
CNumerics::ComputeMeanRateOfStrainMatrix(3, StrainMat, Grad_Vel);
225+
nodes->SetVortex_Tilting(iPoint, StrainMat, Vorticity, Laminar_Viscosity);
221226
}
222227
END_SU2_OMP_FOR
223228
}
@@ -1092,6 +1097,9 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C
10921097
su2double DES_lengthScale = 0.0;
10931098

10941099
switch(kind_hybridRANSLES){
1100+
/*--- Every model is taken from "Development of DDES and IDDES Formulations for the k-ω Shear Stress Transport Model"
1101+
Mikhail S. Gritskevich et al. (DOI:10.1007/s10494-011-9378-4)
1102+
---*/
10951103
case SST_DDES: {
10961104

10971105
const su2double r_d = (eddyVisc + lamVisc) / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10);
@@ -1172,7 +1180,7 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C
11721180
su2double deltaOmega = -1.0;
11731181
su2double vorticityDir[MAXNDIM] = {};
11741182

1175-
for (auto iDim = 0; iDim < 3; iDim++){
1183+
for (auto iDim = 0; iDim < MAXNDIM; iDim++){
11761184
vorticityDir[iDim] = Vorticity[iDim]/VortMag;
11771185
}
11781186

@@ -1183,13 +1191,13 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C
11831191
const auto coord_k = geometry->nodes->GetCoord(kPoint);
11841192

11851193
su2double delta[MAXNDIM] = {};
1186-
// This should only be performed on 3D cases anyway
1187-
for (auto iDim = 0u; iDim < 3; iDim++){
1188-
delta[iDim] = (coord_j[iDim] - coord_k[iDim])/2.0; // Should I divide by 2 as I am interested in the dual volume?
1194+
for (auto iDim = 0u; iDim < MAXNDIM; iDim++){
1195+
// TODO: Should I divide by 2 as I am interested in the dual volume (the edge is split at midpoint)?
1196+
delta[iDim] = (coord_j[iDim] - coord_k[iDim])/2.0;
11891197
}
1190-
su2double l_n_minus_m[3];
1198+
su2double l_n_minus_m[MAXNDIM];
11911199
GeometryToolbox::CrossProduct(delta, vorticityDir, l_n_minus_m);
1192-
deltaOmega = max(deltaOmega, GeometryToolbox::Norm(nDim, l_n_minus_m));
1200+
deltaOmega = max(deltaOmega, GeometryToolbox::Norm(3, l_n_minus_m));
11931201
}
11941202

11951203
// Add to VTM(iPoint) to perform the average

‎SU2_CFD/src/variables/CTurbVariable.cpp‎

Lines changed: 14 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727

2828

2929
#include "../../include/variables/CTurbVariable.hpp"
30+
#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
3031

3132

3233
CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config)
@@ -43,43 +44,29 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned
4344
}
4445

4546

46-
void CTurbVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
47-
const su2double* Vorticity, su2double LaminarViscosity) {
47+
void CTurbVariable::SetVortex_Tilting(unsigned long iPoint, const su2double Strain[3][3],
48+
const su2double* Vorticity, su2double LaminarViscosity) {
4849

49-
su2double Strain[3][3] = {{0,0,0}, {0,0,0}, {0,0,0}}, Omega, StrainDotVort[3], numVecVort[3];
50+
su2double Omega, StrainDotVort[3], numVecVort[3];
5051
su2double numerator, trace0, trace1, denominator;
5152

5253
AD::StartPreacc();
53-
AD::SetPreaccIn(PrimGrad_Flow, nDim+1, nDim);
54+
AD::SetPreaccIn(Strain, nDim, nDim);
5455
AD::SetPreaccIn(Vorticity, 3);
5556
/*--- Eddy viscosity ---*/
5657
AD::SetPreaccIn(muT(iPoint));
5758
/*--- Laminar viscosity --- */
5859
AD::SetPreaccIn(LaminarViscosity);
5960

60-
Strain[0][0] = PrimGrad_Flow[1][0];
61-
Strain[1][0] = 0.5*(PrimGrad_Flow[2][0] + PrimGrad_Flow[1][1]);
62-
Strain[0][1] = 0.5*(PrimGrad_Flow[1][1] + PrimGrad_Flow[2][0]);
63-
Strain[1][1] = PrimGrad_Flow[2][1];
64-
if (nDim == 3){
65-
Strain[0][2] = 0.5*(PrimGrad_Flow[3][0] + PrimGrad_Flow[1][2]);
66-
Strain[1][2] = 0.5*(PrimGrad_Flow[3][1] + PrimGrad_Flow[2][2]);
67-
Strain[2][0] = 0.5*(PrimGrad_Flow[1][2] + PrimGrad_Flow[3][0]);
68-
Strain[2][1] = 0.5*(PrimGrad_Flow[2][2] + PrimGrad_Flow[3][1]);
69-
Strain[2][2] = PrimGrad_Flow[3][2];
70-
}
71-
72-
Omega = sqrt(Vorticity[0]*Vorticity[0] + Vorticity[1]*Vorticity[1]+ Vorticity[2]*Vorticity[2]);
73-
74-
StrainDotVort[0] = Strain[0][0]*Vorticity[0]+Strain[0][1]*Vorticity[1]+Strain[0][2]*Vorticity[2];
75-
StrainDotVort[1] = Strain[1][0]*Vorticity[0]+Strain[1][1]*Vorticity[1]+Strain[1][2]*Vorticity[2];
76-
StrainDotVort[2] = Strain[2][0]*Vorticity[0]+Strain[2][1]*Vorticity[1]+Strain[2][2]*Vorticity[2];
77-
78-
numVecVort[0] = StrainDotVort[1]*Vorticity[2] - StrainDotVort[2]*Vorticity[1];
79-
numVecVort[1] = StrainDotVort[2]*Vorticity[0] - StrainDotVort[0]*Vorticity[2];
80-
numVecVort[2] = StrainDotVort[0]*Vorticity[1] - StrainDotVort[1]*Vorticity[0];
81-
82-
numerator = sqrt(6.0) * sqrt(numVecVort[0]*numVecVort[0] + numVecVort[1]*numVecVort[1] + numVecVort[2]*numVecVort[2]);
61+
Omega = GeometryToolbox::Norm(3, Vorticity);
62+
63+
StrainDotVort[0] = GeometryToolbox::DotProduct(3, Strain[0], Vorticity);
64+
StrainDotVort[1] = GeometryToolbox::DotProduct(3, Strain[1], Vorticity);
65+
StrainDotVort[2] = GeometryToolbox::DotProduct(3, Strain[2], Vorticity);
66+
67+
GeometryToolbox::CrossProduct(StrainDotVort, Vorticity, numVecVort);
68+
69+
numerator = sqrt(6.0) * GeometryToolbox::Norm(3, numVecVort);
8370
trace0 = 3.0*(pow(Strain[0][0],2.0) + pow(Strain[1][1],2.0) + pow(Strain[2][2],2.0));
8471
trace1 = pow(Strain[0][0] + Strain[1][1] + Strain[2][2],2.0);
8572
denominator = pow(Omega, 2.0) * sqrt(trace0-trace1);

0 commit comments

Comments
 (0)