DDES, IDDES and shear-layer-adapted DDES for SST, and fixes to SA EDDES

.github/workflows/regression.yml

@@ -235,7 +235,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c feature_SST_IDDES -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:
@@ -284,7 +284,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests"
+          args: -b ${{github.ref}} -t develop -c feature_SST_IDDES -s ${{matrix.testscript}} -a "--tapetests"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:
@@ -334,7 +334,7 @@ jobs:
           PMIX_MCA_gds: hash
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan"
+          args: -b ${{github.ref}} -t develop -c feature_SST_IDDES -s ${{matrix.testscript}} -a "--tsan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:261011-1107
         with:
@@ -379,7 +379,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:261011-1107
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan"
+          args: -b ${{github.ref}} -t develop -c feature_SST_IDDES -s ${{matrix.testscript}} -a "--asan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:261011-1107
         with:

Common/include/CConfig.hpp

@@ -9836,6 +9836,11 @@ class CConfig {
    */
   unsigned short GetKind_HybridRANSLES(void) const { return Kind_HybridRANSLES; }
 
+  bool IsHybridSST() const {
+    return Kind_HybridRANSLES == SST_DDES || Kind_HybridRANSLES == SST_IDDES ||
+           Kind_HybridRANSLES == SST_SIDDES || Kind_HybridRANSLES == SST_EDDES;
+  }
+
   /*!
    * \brief Get if the LES mode must be enforced.
    * \return TRUE if LES is enforced.

Common/include/geometry/dual_grid/CPoint.hpp

@@ -99,6 +99,7 @@ class CPoint {
   su2vector<unsigned short> nNeighbor; /*!< \brief Number of neighbors, needed by some numerical methods. */
 
   /*--- Closest element on a viscous wall, and distance to it. ---*/
+  su2activematrix WallNormal;                 /*!< \brief Unit normal of the nearest wall element. */
   su2activevector Wall_Distance;              /*!< \brief Distance to the nearest wall. */
   su2vector<int> ClosestWall_Rank;            /*!< \brief Rank of process holding the closest wall element. */
   su2vector<unsigned short> ClosestWall_Zone; /*!< \brief Zone index of closest wall element. */
@@ -887,6 +888,23 @@ class CPoint {
     }
   }
 
+  /*! \brief Gather the nearest wall-element normal for each point. */
+  template <typename NormalsType>
+  void SetWallNormals(const NormalsType& normals) {
+    WallNormal.resize(GlobalIndex.size(), nDim) = su2double(0.0);
+    for (unsigned long iPoint = 0; iPoint < GlobalIndex.size(); ++iPoint) {
+      const auto rank = ClosestWall_Rank[iPoint];
+      if (rank < 0) continue;
+      const auto& normal =
+          normals[rank][ClosestWall_Zone[iPoint]][ClosestWall_Marker[iPoint]][ClosestWall_Elem[iPoint]];
+      for (auto iDim = 0u; iDim < nDim; ++iDim) WallNormal(iPoint, iDim) = normal[iDim];
+    }
+  }
+
+  inline const su2double* GetWallNormal(unsigned long iPoint) const {
+    return WallNormal.empty() ? nullptr : WallNormal[iPoint];
+  }
+
   /*!
    * \brief Set wall roughnesses according to stored closest wall information.
    * \param[in] roughness - Mapping [rank][zone][marker] -> roughness

Common/include/option_structure.hpp

@@ -1196,7 +1196,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   const bool sst_compWilcox = IsPresent(SST_OPTIONS::COMP_Wilcox);
   const bool sst_compSarkar = IsPresent(SST_OPTIONS::COMP_Sarkar);
   const bool sst_dll = IsPresent(SST_OPTIONS::DLL);
-
   if (sst_1994 && sst_2003) {
     SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
   } else if (sst_2003) {
@@ -1702,17 +1701,25 @@ static const MapType<std::string, WINDOW_FUNCTION> Window_Map = {
  */
 enum ENUM_HYBRIDRANSLES {
   NO_HYBRIDRANSLES = 0,  /*!< \brief No turbulence model. */
-  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_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). */
+  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). */
+  SST_EDDES        = 8   /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
 };
 static const MapType<std::string, ENUM_HYBRIDRANSLES> HybridRANSLES_Map = {
   MakePair("NONE", NO_HYBRIDRANSLES)
   MakePair("SA_DES", SA_DES)
   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)
+  MakePair("SST_EDDES", SST_EDDES)
 };
 
 /*!

Common/include/toolboxes/geometry_toolbox.hpp

@@ -60,6 +60,51 @@ inline void PointPointReflect(Int nDim, const T* a, const T* b, T* d) {
   for (Int i = 0; i < nDim; i++) d[i] = 2 * b[i] - a[i];
 }
 
+/*!
+ * \brief Compute the mean rate of strain matrix.
+ * \details The parameter primvargrad can be e.g. PrimVar_Grad_i or Mean_GradPrimVar.
+ * \param[in] nDim - 2 or 3
+ * \param[out] rateofstrain - Rate of strain matrix
+ * \param[in] velgrad - A velocity gradient matrix.
+ * \tparam Mat1 - any type that supports the [][] interface
+ * \tparam Mat2 - any type that supports the [][] interface
+ */
+template <class Mat1, class Mat2>
+inline void ComputeMeanRateOfStrainMatrix(size_t nDim, Mat1& rateofstrain, const Mat2& velgrad) {
+  /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */
+
+  if (nDim == 3) {
+    rateofstrain[0][0] = velgrad[0][0];
+    rateofstrain[1][1] = velgrad[1][1];
+    rateofstrain[2][2] = velgrad[2][2];
+    rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]);
+    rateofstrain[0][2] = 0.5 * (velgrad[0][2] + velgrad[2][0]);
+    rateofstrain[1][2] = 0.5 * (velgrad[1][2] + velgrad[2][1]);
+    rateofstrain[1][0] = rateofstrain[0][1];
+    rateofstrain[2][1] = rateofstrain[1][2];
+    rateofstrain[2][0] = rateofstrain[0][2];
+  } else {  // nDim==2
+    rateofstrain[0][0] = velgrad[0][0];
+    rateofstrain[1][1] = velgrad[1][1];
+    rateofstrain[2][2] = 0.0;
+    rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]);
+    rateofstrain[0][2] = 0.0;
+    rateofstrain[1][2] = 0.0;
+    rateofstrain[1][0] = rateofstrain[0][1];
+    rateofstrain[2][1] = rateofstrain[1][2];
+    rateofstrain[2][0] = rateofstrain[0][2];
+  }
+}
+
+/*! \brief Magnitude of an edge projected onto a unit wall normal. */
+template <class T, typename Int>
+inline T WallNormalEdgeLength(Int nDim, const T* a, const T* b, const T* normal) {
+  T projection(0);
+  for (Int i = 0; i < nDim; ++i) projection += (b[i] - a[i]) * normal[i];
+  using std::abs;
+  return abs(projection);
+}
+
 /*! \return a.b */
 template <class T, typename Int>
 inline T DotProduct(Int nDim, const T* a, const T* b) {

Common/src/CConfig.cpp

@@ -3124,7 +3124,7 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Window (weight) function for the cost-functional in the reverse sweep */
   addEnumOption("WINDOW_FUNCTION", Kind_WindowFct, Window_Map, WINDOW_FUNCTION::SQUARE);
 
-  /* DESCRIPTION: DES Constant */
+  /* DESCRIPTION: DES Constant of the SA hybrid models. The SST hybrid models blend their own constants with F1. */
   addDoubleOption("DES_CONST", Const_DES, 0.65);
 
   /* DESCRIPTION: SBS lengthscale coefficient */
@@ -3145,7 +3145,7 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Specify if the Stochastic Backscatter Model must be activated */
   addBoolOption("STOCHASTIC_BACKSCATTER", SBSParam.StochasticBackscatter, false);
 
-  /* DESCRIPTION: Specify if the LES mode must be enforced */
+  /* DESCRIPTION: Specify if the LES mode must be enforced (SA hybrid models only) */
   addBoolOption("ENFORCE_LES", enforceLES, false);
 
   /* DESCRIPTION: Specify if the stochastic source term must be included in the turbulence model equation */
@@ -3166,7 +3166,7 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Shielding function lower threshold for application of Stochastic Backscatter Model */
   addDoubleOption("SBS_FD_LOWER_THRESHOLD", SBSParam.stochFdThreshold, 0.9);
 
-  /* DESCRIPTION: Filter width for LES (if negative, it is computed based on the local cell size) */
+  /* DESCRIPTION: Filter width for LES (if negative, it is computed based on the local cell size, SA hybrid models only) */
   addDoubleOption("LES_FILTER_WIDTH", LES_FilterWidth, -1.0);
 
   /* DESCRIPTION: Roe with low dissipation for unsteady flows */
@@ -4149,6 +4149,52 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
                    "TIME_STEPPING, DUAL_TIME_STEPPING-1ST_ORDER or DUAL_TIME_STEPPING-2ND_ORDER", CURRENT_FUNCTION);
   }
 
+  if (Kind_HybridRANSLES != NO_HYBRIDRANSLES && !Time_Domain) {
+    SU2_MPI::Error("Hybrid RANS/LES models (HYBRID_RANSLES) require an unsteady simulation, set TIME_DOMAIN= YES.",
+                   CURRENT_FUNCTION);
+  }
+
+  /*--- The SA_* hybrid models are built on the SA model and the SST_* ones on the SST model. ---*/
+  const bool hybridSA = Kind_HybridRANSLES == SA_DES || Kind_HybridRANSLES == SA_DDES ||
+                        Kind_HybridRANSLES == SA_ZDES || Kind_HybridRANSLES == SA_EDDES;
+  const bool hybridSST = IsHybridSST();
+  if ((hybridSA && Kind_Turb_Model != TURB_MODEL::SA) || (hybridSST && Kind_Turb_Model != TURB_MODEL::SST)) {
+    SU2_MPI::Error(string("The hybrid RANS/LES model (HYBRID_RANSLES) must match the turbulence model (KIND_TURB_MODEL): ") +
+                   "SA_* options require SA, SST_* options require SST.", CURRENT_FUNCTION);
+  }
+
+  /*--- The stochastic backscatter model adds the Langevin equations to the SA model (turbulence variables 1 to 3). ---*/
+  if (SBSParam.StochasticBackscatter && hybridSST) {
+    SU2_MPI::Error(string("The stochastic backscatter model (STOCHASTIC_BACKSCATTER= YES) is only available with the SA ") +
+                   "hybrid RANS/LES models (HYBRID_RANSLES= SA_DES, SA_DDES, SA_ZDES or SA_EDDES).", CURRENT_FUNCTION);
+  }
+
+  /*--- How the SST compressibility corrections should act in the LES region of the hybrid models is not validated yet
+   (e.g. Forsythe et al. 2002 keep them only in the RANS length scale). Stop rather than give unvalidated results. ---*/
+  if (hybridSST && (sstParsedOptions.compWilcox || sstParsedOptions.compSarkar)) {
+    SU2_MPI::Error(string("The SST compressibility corrections (COMPRESSIBILITY-WILCOX, COMPRESSIBILITY-SARKAR) are not ") +
+                   "supported yet with the SST hybrid RANS/LES models (HYBRID_RANSLES= SST_*).", CURRENT_FUNCTION);
+  }
+
+  /*--- The SST hybrid models compute their length scale from the grid with the constant C_DES = 0.78 F1 + 0.61 (1 - F1)
+   (Gritskevich et al. 2012), they do not use these options of the SA hybrid models. ---*/
+  if (hybridSST) {
+    if (enforceLES) {
+      SU2_MPI::Error("ENFORCE_LES= YES is only available with the SA hybrid RANS/LES models (HYBRID_RANSLES= SA_*).",
+                     CURRENT_FUNCTION);
+    }
+    if (LES_FilterWidth > 0.0) {
+      SU2_MPI::Error(string("A user-specified LES_FILTER_WIDTH is only available with the SA hybrid RANS/LES models ") +
+                     "(HYBRID_RANSLES= SA_*). The SST hybrid models compute the filter width from the grid.", CURRENT_FUNCTION);
+    }
+    if (OptionIsSet("DES_CONST")) {
+      if (rank == MASTER_NODE) {
+        cout << "WARNING: DES_CONST is ignored by SST hybrid RANS/LES models, which use "
+                "C_DES = 0.78 F1 + 0.61 (1 - F1)." << endl;
+      }
+    }
+  }
+
   if (Time_Domain){
     Delta_UnstTime = Time_Step;
 
@@ -6847,13 +6893,21 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         cout << "Hybrid RANS/LES: ";
         switch (Kind_HybridRANSLES) {
           case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break;
-          case SA_DES:   cout << "Detached Eddy Simulation (DES97) " << endl; break;
-          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 SA_DES:         cout << "Detached Eddy Simulation (DES97) " << endl; break;
+          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)" << 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;
+          case SST_EDDES:      cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
         }
         if (Kind_HybridRANSLES != NO_HYBRIDRANSLES) {
           if (LES_FilterWidth > 0.0) cout << "User-specified LES filter width: " << LES_FilterWidth << endl;
+          if (IsHybridSST())
+            cout << "DES constant blended with F1: C_DES = 0.78 F1 + 0.61 (1 - F1)." << endl;
+          else
+            cout << "DES constant: " << Const_DES << endl;
           cout << "Stochastic Backscatter: ";
           if (SBSParam.StochasticBackscatter) {
             cout << "ON" << endl;

Common/src/geometry/CGeometry.cpp

@@ -4542,6 +4542,52 @@ su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, co
   const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint);
   return 2 * Vol / GeometryToolbox::Norm(3, Normal);
 }
+/*! \brief Average the vertex normals of a wall element and normalize the result. */
+std::array<su2double, 3> WallElementNormal(const CGeometry* geometry, unsigned short iMarker, unsigned long iElem) {
+  std::array<su2double, 3> normal = {};
+  const auto* element = geometry->bound[iMarker][iElem];
+  for (auto iNode = 0u; iNode < element->GetnNodes(); ++iNode) {
+    const auto iPoint = element->GetNode(iNode);
+    const auto iVertex = geometry->nodes->GetVertex(iPoint, iMarker);
+    for (auto iDim = 0u; iDim < geometry->GetnDim(); ++iDim)
+      normal[iDim] += geometry->vertex[iMarker][iVertex]->GetNormal(iDim);
+  }
+  for (auto& component : normal) component /= element->GetnNodes();
+  const auto magnitude = GeometryToolbox::Norm(3, normal.data());
+  if (magnitude > 0.0)
+    for (auto& component : normal) component /= magnitude;
+  return normal;
+}
+
+/*! \brief Collect unit wall-element normals for each marker of one zone. */
+su2vector<su2matrix<su2double>> CollectWallNormals(const CGeometry* geometry, const CConfig* config) {
+  su2vector<su2matrix<su2double>> normals;
+  normals.resize(geometry->GetnMarker());
+  for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); ++iMarker) {
+    if (!config->GetViscous_Wall(iMarker)) {
+      normals[iMarker].resize(1, 3) = su2double(0.0);
+      continue;
+    }
+    normals[iMarker].resize(geometry->GetnElem_Bound(iMarker), 3);
+    for (auto iElem = 0ul; iElem < geometry->GetnElem_Bound(iMarker); ++iElem) {
+      const auto normal = WallElementNormal(geometry, iMarker, iElem);
+      for (auto iDim = 0u; iDim < 3; ++iDim) normals[iMarker](iElem, iDim) = normal[iDim];
+    }
+  }
+  return normals;
+}
+
+/*! \brief Adapt one marker's normal matrix to the recursive NdFlattener interface. */
+auto WallNormalMatrix(const su2matrix<su2double>& normals) {
+  return make_pair(normals.rows(), [&normals](auto iElem) {
+    return make_pair(normals.cols(), [&normals, iElem](auto iDim) { return normals(iElem, iDim); });
+  });
+}
+
+/*! \brief Adapt a zone's marker normals without copying the matrices. */
+auto WallNormalMarkers(const su2vector<su2matrix<su2double>>& normals) {
+  return make_pair(normals.size(), [&normals](auto iMarker) { return WallNormalMatrix(normals[iMarker]); });
+}
 }  // namespace
 
 void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container,
@@ -4612,6 +4658,27 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo
       }
     }
 
+    /*--- IDDES needs the normal of the nearest wall, including walls on other ranks and zones. ---*/
+    bool needWallNormals = false;
+    for (int iZone = 0; iZone < nZone; ++iZone) {
+      const auto hybrid = config_container[iZone]->GetKind_HybridRANSLES();
+      needWallNormals = needWallNormals || hybrid == SST_IDDES || hybrid == SST_SIDDES;
+    }
+    if (needWallNormals) {
+      su2vector<su2vector<su2matrix<su2double>>> wallNormals;
+      wallNormals.resize(nZone);
+      for (int iZone = 0; iZone < nZone; ++iZone)
+        wallNormals[iZone] = CollectWallNormals(geometry_container[iZone][iInst][MESH_0], config_container[iZone]);
+      auto normal_i = make_pair(nZone, [&wallNormals](auto iZone) { return WallNormalMarkers(wallNormals[iZone]); });
+      NdFlattener<4> localNormals(normal_i);
+      NdFlattener<5> globalNormals(Nd_MPI_Environment(), localNormals);
+      for (int iZone = 0; iZone < nZone; ++iZone) {
+        const auto hybrid = config_container[iZone]->GetKind_HybridRANSLES();
+        if (hybrid == SST_IDDES || hybrid == SST_SIDDES)
+          geometry_container[iZone][iInst][MESH_0]->nodes->SetWallNormals(globalNormals);
+      }
+    }
+
     for (int iZone = 0; iZone < nZone; iZone++) {
       /*--- When recording for a specific zone, only compute nearest-neighbor distances (which read vertex
        *    normals) for that zone. Reading normals from other zones at this tape position would create