Implementation of Simplified LM Transition model

.gitmodules

@@ -1,6 +1,6 @@
 [submodule "externals/codi"]
 	path = externals/codi
-	url = https://github.com/scicompkl/CoDiPack
+	url = https://github.com/SciCompKL/CoDiPack
 	ignore = dirty
 [submodule "externals/medi"]
 	path = externals/medi

Common/include/geometry/dual_grid/CPoint.hpp

@@ -111,6 +111,7 @@ class CPoint {
   su2activevector Curvature;          /*!< \brief Value of the surface curvature (SU2_GEO). */
   su2activevector MaxLength;          /*!< \brief The maximum cell-center to cell-center length. */
   su2activevector RoughnessHeight;    /*!< \brief Roughness of the nearest wall. */
+  su2activematrix Normals;            /*!< \brief Normal of the nearest wall element. */
 
   su2matrix<AD::Identifier>
       AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector before solver iteration. */
@@ -484,6 +485,27 @@ class CPoint {
   }
   inline void SetWall_Distance(unsigned long iPoint, su2double distance) { Wall_Distance(iPoint) = distance; }
 
+  /*!
+   * \brief Get the index of the closest wall element.
+   * \param[in] iPoint - Index of the point.
+   * \param[out] ClosestWall_Elem - ID of the closest element on a wall boundary.
+   */
+  inline unsigned long GetClosestWall_Elem(unsigned long iPoint) { return ClosestWall_Elem(iPoint); }
+
+  /*!
+   * \brief Get the marker of the closest wall marker.
+   * \param[in] iPoint - Index of the point.
+   * \param[out] ClosestWall_Marker - MarkerID of the closest wall boundary.
+   */
+  inline unsigned long GetClosestWall_Marker(unsigned long iPoint) { return ClosestWall_Marker(iPoint); }
+
+  /*!
+   * \brief Get the rank of the closest wall marker.
+   * \param[in] iPoint - Index of the point.
+   * \param[out] ClosestWall_Rank - RankID of the closest wall boundary.
+   */
+  inline unsigned long GetClosestWall_Rank(unsigned long iPoint) { return ClosestWall_Rank(iPoint); }
+
   /*!
    * \brief Get the value of the distance to the nearest wall.
    * \param[in] iPoint - Index of the point.
@@ -506,6 +528,23 @@ class CPoint {
    */
   inline su2double GetRoughnessHeight(unsigned long iPoint) const { return RoughnessHeight(iPoint); }
 
+  /*!
+   * \brief Set the value of the normal of the nearest wall element.
+   * \param[in] iPoint - Index of the point.
+   * \param[in] normal - Value of the normal.
+   */
+  template <typename Normals_type>
+  inline void SetNormal(unsigned long iPoint, Normals_type const& normal) {
+    for (unsigned long iDim = 0; iDim < nDim; iDim++) Normals(iPoint, iDim) = normal[iDim];
+  }
+
+  /*!
+   * \brief Set the value of the normal of the nearest wall element.
+   * \param[in] iPoint - Index of the point.
+   * \return normal to the normal of the nearest wall element.
+   */
+  inline su2double* GetNormal(unsigned long iPoint) { return Normals[iPoint]; }
+
   /*!
    * \brief Set the value of the distance to a sharp edge.
    * \param[in] iPoint - Index of the point.
@@ -885,4 +924,19 @@ class CPoint {
       }
     }
   }
+
+  /*!
+   * \brief Set wall normal according to stored closest wall information.
+   * \param[in] normals - Mapping [rank][zone][marker][element] -> normal
+   */
+  template <typename Normals_type>
+  void SetWallNormals(Normals_type const& normals) {
+    for (unsigned long iPoint = 0; iPoint < GlobalIndex.size(); ++iPoint) {
+      auto rankID = ClosestWall_Rank[iPoint];
+      auto zoneID = ClosestWall_Zone[iPoint];
+      auto markerID = ClosestWall_Marker[iPoint];
+      auto elementID = ClosestWall_Elem[iPoint];
+      if (rankID >= 0) SetNormal(iPoint, normals[rankID][zoneID][markerID][elementID]);
+    }
+  }
 };

Common/include/option_structure.hpp

@@ -1216,27 +1216,37 @@ static const MapType<std::string, TURB_TRANS_MODEL> Trans_Model_Map = {
  * \brief LM Options
  */
 enum class LM_OPTIONS {
-  NONE,         /*!< \brief No option / default. */
-  LM2015,       /*!< \brief Cross-flow corrections. */
-  MALAN,        /*!< \brief Kind of transition correlation model (Malan). */
-  SULUKSNA,     /*!< \brief Kind of transition correlation model (Suluksna). */
-  KRAUSE,       /*!< \brief Kind of transition correlation model (Krause). */
-  KRAUSE_HYPER, /*!< \brief Kind of transition correlation model (Krause hypersonic). */
-  MEDIDA_BAEDER,/*!< \brief Kind of transition correlation model (Medida-Baeder). */
-  MEDIDA,       /*!< \brief Kind of transition correlation model (Medida). */
-  MENTER_LANGTRY,   /*!< \brief Kind of transition correlation model (Menter-Langtry). */
-  DEFAULT       /*!< \brief Kind of transition correlation model (Menter-Langtry if SST, MALAN if SA). */
+  NONE,           /*!< \brief No option / default. */
+  CROSSFLOW,      /*!< \brief Cross-flow corrections. */
+  SLM,            /*!< \brief Simplified version. */
+  PRODLIM,        /*!< \brief Add production term to Pk. */
+  MALAN,          /*!< \brief Kind of transition correlation model (Malan). */
+  SULUKSNA,       /*!< \brief Kind of transition correlation model (Suluksna). */
+  KRAUSE,         /*!< \brief Kind of transition correlation model (Krause). */
+  KRAUSE_HYPER,   /*!< \brief Kind of transition correlation model (Krause hypersonic). */
+  MEDIDA_BAEDER,  /*!< \brief Kind of transition correlation model (Medida-Baeder). */
+  MEDIDA,         /*!< \brief Kind of transition correlation model (Medida). */
+  MENTER_LANGTRY, /*!< \brief Kind of transition correlation model (Menter-Langtry). */
+  MENTER_SLM,     /*!< \brief Kind of transition correlation model (Menter Simplified LM model). */
+  CODER_SLM,      /*!< \brief Kind of transition correlation model (Coder Simplified LM model). */
+  MOD_EPPLER_SLM, /*!< \brief Kind of transition correlation model (Modified Eppler Simplified LM model). */
+  DEFAULT         /*!< \brief Kind of transition correlation model (Menter-Langtry if SST, MALAN if SA). */
 };
 
 static const MapType<std::string, LM_OPTIONS> LM_Options_Map = {
   MakePair("NONE", LM_OPTIONS::NONE)
-  MakePair("LM2015", LM_OPTIONS::LM2015)
+  MakePair("CROSSFLOW", LM_OPTIONS::CROSSFLOW)
+  MakePair("SLM", LM_OPTIONS::SLM)
+  MakePair("PRODLIM", LM_OPTIONS::PRODLIM)
   MakePair("MALAN", LM_OPTIONS::MALAN)
   MakePair("SULUKSNA", LM_OPTIONS::SULUKSNA)
   MakePair("KRAUSE", LM_OPTIONS::KRAUSE)
   MakePair("KRAUSE_HYPER", LM_OPTIONS::KRAUSE_HYPER)
   MakePair("MEDIDA_BAEDER", LM_OPTIONS::MEDIDA_BAEDER)
   MakePair("MENTER_LANGTRY", LM_OPTIONS::MENTER_LANGTRY)
+  MakePair("MENTER_SLM", LM_OPTIONS::MENTER_SLM)
+  MakePair("CODER_SLM", LM_OPTIONS::CODER_SLM)
+  MakePair("MOD_EPPLER_SLM", LM_OPTIONS::MOD_EPPLER_SLM)
   MakePair("DEFAULT", LM_OPTIONS::DEFAULT)
 };
 
@@ -1254,13 +1264,26 @@ enum class TURB_TRANS_CORRELATION {
   DEFAULT       /*!< \brief Kind of transition correlation model (Menter-Langtry if SST, MALAN if SA). */
 };
 
+/*!
+ * \brief Types of transition correlations for Simplified LM model
+ */
+enum class TURB_TRANS_CORRELATION_SLM {
+  MENTER_SLM,        /*!< \brief Kind of transition correlation model (Menter Simplified LM model). */
+  CODER_SLM,        /*!< \brief Kind of transition correlation model (Coder Simplified LM model). */
+  MOD_EPPLER_SLM,        /*!< \brief Kind of transition correlation model (Modified Eppler Simplified LM model). */
+  DEFAULT       /*!< \brief Kind of transition correlation model. */
+};
+
 /*!
  * \brief Structure containing parsed LM options.
  */
 struct LM_ParsedOptions {
   LM_OPTIONS version = LM_OPTIONS::NONE;  /*!< \brief LM base model. */
-  bool LM2015 = false;                    /*!< \brief Use cross-flow corrections. */
+  bool SLM = false;                       /*!< \brief Use simplified version. */
+  bool ProdLim = false;                   /*!< \brief Add production term to Pk. */
+  bool CrossFlow = false;                    /*!< \brief Use cross-flow corrections. */
   TURB_TRANS_CORRELATION Correlation = TURB_TRANS_CORRELATION::DEFAULT;
+  TURB_TRANS_CORRELATION_SLM Correlation_SLM = TURB_TRANS_CORRELATION_SLM::DEFAULT;
 };
 
 /*!
@@ -1278,7 +1301,10 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh
     return std::find(LM_Options, lm_options_end, option) != lm_options_end;
   };
 
-  LMParsedOptions.LM2015 = IsPresent(LM_OPTIONS::LM2015);
+  LMParsedOptions.SLM = IsPresent(LM_OPTIONS::SLM);
+  LMParsedOptions.ProdLim = IsPresent(LM_OPTIONS::PRODLIM);
+
+  LMParsedOptions.CrossFlow = IsPresent(LM_OPTIONS::CROSSFLOW);
 
   int NFoundCorrelations = 0;
   if (IsPresent(LM_OPTIONS::MALAN)) {
@@ -1310,10 +1336,28 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh
     NFoundCorrelations++;
   }
 
+  int NFoundCorrelations_SLM = 0;
+  if (IsPresent(LM_OPTIONS::MENTER_SLM)) {
+    LMParsedOptions.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::MENTER_SLM;
+    NFoundCorrelations_SLM++;
+  }
+  if (IsPresent(LM_OPTIONS::CODER_SLM)) {
+    LMParsedOptions.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::CODER_SLM;
+    NFoundCorrelations_SLM++;
+  }
+  if (IsPresent(LM_OPTIONS::MOD_EPPLER_SLM)) {
+    LMParsedOptions.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::MOD_EPPLER_SLM;
+    NFoundCorrelations_SLM++;
+  }
+
   if (NFoundCorrelations > 1) {
     SU2_MPI::Error("Two correlations selected for LM_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
   }
 
+  if (NFoundCorrelations_SLM > 1) {
+    SU2_MPI::Error("Two SLM correlations selected for LM_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
+  }
+
   if (LMParsedOptions.Correlation == TURB_TRANS_CORRELATION::DEFAULT){
     if (Kind_Turb_Model == TURB_MODEL::SST) {
       LMParsedOptions.Correlation = TURB_TRANS_CORRELATION::MENTER_LANGTRY;
@@ -1322,6 +1366,10 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh
     }
   }
 
+  if (LMParsedOptions.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::DEFAULT){
+    LMParsedOptions.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::MENTER_SLM;
+  }
+
   return LMParsedOptions;
 }
 

Common/src/CConfig.cpp

@@ -3583,9 +3583,9 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) {
     lmParsedOptions = ParseLMOptions(LM_Options, nLM_Options, rank, Kind_Turb_Model);
 
-    /*--- Check if problem is 2D and LM2015 has been selected ---*/
-    if (lmParsedOptions.LM2015 && val_nDim == 2) {
-      SU2_MPI::Error("LM2015 is available only for 3D problems", CURRENT_FUNCTION);
+    /*--- Check if problem is 2D and CrossFlow has been selected ---*/
+    if (lmParsedOptions.CrossFlow && val_nDim == 2) {
+      SU2_MPI::Error("Cross-flow corrections are available only for 3D problems", CURRENT_FUNCTION);
     }
   }
 
@@ -6410,33 +6410,62 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         switch (Kind_Trans_Model) {
           case TURB_TRANS_MODEL::NONE:  break;
           case TURB_TRANS_MODEL::LM: {
-            cout << "Transition model: Langtry and Menter's 4 equation model";
-            if (lmParsedOptions.LM2015) {
-              cout << " w/ cross-flow corrections (2015)" << endl;
+            int NTurbEqs = 0;
+            switch (Kind_Turb_Model) {
+              case TURB_MODEL::SA: NTurbEqs = 1;  break;
+              case TURB_MODEL::SST: NTurbEqs = 2;  break;
+              case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but LM transition model is active.", CURRENT_FUNCTION); break;
+            }
+            if (!lmParsedOptions.SLM) {
+              int NEquations = 2;
+              cout << "Transition model: Langtry and Menter's "<< NEquations+NTurbEqs <<" equation model";
             } else {
-              cout << " (2009)" << endl;
+              int NEquations = 1;
+              cout << "Transition model: Simplified Langtry and Menter's "<< NEquations+NTurbEqs <<" equation model";
+            }
+            if (lmParsedOptions.CrossFlow) {
+              cout << " w/ cross-flow corrections";
+              if (!lmParsedOptions.SLM) {
+                cout << " (2015)";
+              }
+              cout << endl;
+            } else {
+              if (!lmParsedOptions.SLM) {
+                cout << " (2009)" << endl;
+              }
             }
             break;
           }
         }
-        if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) {
 
-          cout << "Correlation Functions: ";
-          switch (lmParsedOptions.Correlation) {
-            case TURB_TRANS_CORRELATION::MALAN: cout << "Malan et al. (2009)" << endl;  break;
-            case TURB_TRANS_CORRELATION::SULUKSNA: cout << "Suluksna et al. (2009)" << endl;  break;
-            case TURB_TRANS_CORRELATION::KRAUSE: cout << "Krause et al. (2008)" << endl;  break;
-            case TURB_TRANS_CORRELATION::KRAUSE_HYPER: cout << "Krause et al. (2008, paper)" << endl;  break;
-            case TURB_TRANS_CORRELATION::MEDIDA_BAEDER: cout << "Medida and Baeder (2011)" << endl;  break;
-            case TURB_TRANS_CORRELATION::MEDIDA: cout << "Medida PhD (2014)" << endl;  break;
-            case TURB_TRANS_CORRELATION::MENTER_LANGTRY: cout << "Menter and Langtry (2009)" << endl;  break;
-            case TURB_TRANS_CORRELATION::DEFAULT:
-              switch (Kind_Turb_Model) {
-                case TURB_MODEL::SA: cout << "Malan et al. (2009)" << endl;  break;
-                case TURB_MODEL::SST: cout << "Menter and Langtry (2009)" << endl;  break;
-                case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but LM transition model is active.", CURRENT_FUNCTION); break;
-              }
-              break;
+        if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) {
+          if (!lmParsedOptions.SLM){
+            cout << "Correlation Functions: ";
+            switch (lmParsedOptions.Correlation) {
+              case TURB_TRANS_CORRELATION::MALAN: cout << "Malan et al. (2009)" << endl;  break;
+              case TURB_TRANS_CORRELATION::SULUKSNA: cout << "Suluksna et al. (2009)" << endl;  break;
+              case TURB_TRANS_CORRELATION::KRAUSE: cout << "Krause et al. (2008)" << endl;  break;
+              case TURB_TRANS_CORRELATION::KRAUSE_HYPER: cout << "Krause et al. (2008, paper)" << endl;  break;
+              case TURB_TRANS_CORRELATION::MEDIDA_BAEDER: cout << "Medida and Baeder (2011)" << endl;  break;
+              case TURB_TRANS_CORRELATION::MEDIDA: cout << "Medida PhD (2014)" << endl;  break;
+              case TURB_TRANS_CORRELATION::MENTER_LANGTRY: cout << "Menter and Langtry (2009)" << endl;  break;
+              case TURB_TRANS_CORRELATION::DEFAULT:
+                switch (Kind_Turb_Model) {
+                  case TURB_MODEL::SA: cout << "Malan et al. (2009)" << endl;  break;
+                  case TURB_MODEL::SST: cout << "Menter and Langtry (2009)" << endl;  break;
+                  case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but LM transition model is active.", CURRENT_FUNCTION); break;
+                }
+                break;
+            }
+          }
+          else {
+            cout << "Correlation Functions for Simplified LM model: ";
+            switch (lmParsedOptions.Correlation_SLM) {
+              case TURB_TRANS_CORRELATION_SLM::CODER_SLM: cout << "Coder et al. (2012)" << endl;  break;
+              case TURB_TRANS_CORRELATION_SLM::MOD_EPPLER_SLM: cout << "Modified Eppler (from Coder et al. 2012)" << endl;  break;
+              case TURB_TRANS_CORRELATION_SLM::MENTER_SLM:
+              case TURB_TRANS_CORRELATION_SLM::DEFAULT: cout << "Menter et al. (2015)" << endl;  break;
+            }
           }
         }
         cout << "Hybrid RANS/LES: ";

Common/src/geometry/CGeometry.cpp

@@ -4207,5 +4207,78 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo
         END_SU2_OMP_FOR
       }
     }
+
+    MAIN_SOLVER kindSolver = config_container[ZONE_0]->GetKind_Solver();
+    if (!wallDistanceNeeded[ZONE_0] || kindSolver == MAIN_SOLVER::FEM_LES || kindSolver == MAIN_SOLVER::FEM_RANS) {
+      continue;
+    } else {
+      su2vector<su2vector<su2matrix<su2double>>> WallNormal_container;
+      WallNormal_container.resize(nZone) = su2vector<su2matrix<su2double>>();
+      for (int iZone = 0; iZone < nZone; iZone++) {
+        const CConfig* config = config_container[iZone];
+        const CGeometry* geometry = geometry_container[iZone][iInst][MESH_0];
+        WallNormal_container[iZone].resize(geometry->GetnMarker());
+        for (auto iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) {
+          if (config->GetViscous_Wall(iMarker)) {
+            WallNormal_container[iZone][iMarker].resize(geometry->GetnElem_Bound(iMarker), 3);
+
+            for (auto iElem = 0u; iElem < geometry->GetnElem_Bound(iMarker); iElem++) {
+              su2vector<su2double> NormalHere;
+              NormalHere.resize(3) = su2double(0.0);
+
+              for (unsigned short iNode = 0; iNode < geometry->bound[iMarker][iElem]->GetnNodes(); iNode++) {
+                // Extract global coordinate of the node
+                unsigned long iPointHere = geometry->bound[iMarker][iElem]->GetNode(iNode);
+                long iVertexHere = geometry->nodes->GetVertex(iPointHere, iMarker);
+                for (auto iDim = 0u; iDim < 3; iDim++)
+                  NormalHere[iDim] += geometry->vertex[iMarker][iVertexHere]->GetNormal(iDim);
+              }
+
+              for (auto iDim = 0u; iDim < 3; iDim++) NormalHere[iDim] /= geometry->bound[iMarker][iElem]->GetnNodes();
+
+              su2double NormalMag = 0.0;
+              for (auto iDim = 0u; iDim < 3; iDim++) NormalMag += NormalHere[iDim] * NormalHere[iDim];
+              NormalMag = sqrt(NormalMag);
+
+              for (auto iDim = 0u; iDim < 3; iDim++) NormalHere[iDim] /= NormalMag;
+
+              for (auto iDim = 0u; iDim < 3; iDim++)
+                WallNormal_container[iZone][iMarker](iElem, iDim) = NormalHere[iDim];
+            }
+          } else {
+            WallNormal_container[iZone][iMarker].resize(1, 3) = su2double(0.0);
+          }
+        }
+      }
+
+      auto normal_i =
+          make_pair(nZone, [config_container, geometry_container, iInst, WallNormal_container](unsigned long iZone) {
+            const CConfig* config = config_container[iZone];
+            const CGeometry* geometry = geometry_container[iZone][iInst][MESH_0];
+            const auto nMarker = geometry->GetnMarker();
+            const auto WallNormal = WallNormal_container[iZone];
+
+            return make_pair(nMarker, [config, geometry, WallNormal](unsigned long iMarker) {
+              auto nElem_Bou = geometry->GetnElem_Bound(iMarker);
+              if (!config->GetViscous_Wall(iMarker)) nElem_Bou = 1;
+
+              return make_pair(nElem_Bou, [WallNormal, iMarker](unsigned long iElem) {
+                const auto dimensions = 3;
+
+                return make_pair(dimensions, [WallNormal, iMarker, iElem](unsigned short iDim) {
+                  return WallNormal[iMarker](iElem, iDim);
+                });
+              });
+            });
+          });
+
+      NdFlattener<4> Normals_Local(normal_i);
+      NdFlattener<5> Normals_global(Nd_MPI_Environment(), Normals_Local);
+
+      // use it to update roughnesses
+      for (int jZone = 0; jZone < nZone; jZone++) {
+        geometry_container[jZone][iInst][MESH_0]->nodes->SetWallNormals(Normals_global);
+      }
+    }
   }
 }