[WIP] Implementation of Simplified LM Transition model

Common/include/geometry/dual_grid/CPoint.hpp

@@ -108,10 +108,11 @@ class CPoint {
   su2vector<unsigned long>
       ClosestWall_Elem; /*!< \brief Element index of closest wall element, for givenrank, zone and marker index. */
 
-  su2activevector SharpEdge_Distance; /*!< \brief Distance to a sharp edge. */
-  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. */
+  su2activevector SharpEdge_Distance;     /*!< \brief Distance to a sharp edge. */
+  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<int> AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector. */
   su2matrix<int>
@@ -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,24 @@ 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.
@@ -906,4 +946,22 @@ 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

@@ -1215,27 +1215,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)
 };
 
@@ -1253,13 +1263,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;
 };
 
 /*!
@@ -1277,7 +1300,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)) {
@@ -1309,6 +1335,20 @@ 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);
   }
@@ -1321,6 +1361,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

@@ -3495,9 +3495,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);
     }
   }
 
@@ -6230,11 +6230,27 @@ 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" << endl;
+            } else {
+              if (!lmParsedOptions.SLM) {
+                cout << " (2009)" << endl;
+              } else {
+                cout << " (2015)" << endl;
+              }
             }
             break;
           }
@@ -6258,6 +6274,13 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
               }
               break;
           }
+          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: ";
         switch (Kind_HybridRANSLES) {

Common/src/geometry/CGeometry.cpp

@@ -3969,5 +3969,81 @@
         }
       }
     }
+
+
+    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);
+
+          // cout << "geometry->nVertex[iMarker] = " << geometry->nVertex[iMarker] << endl;
+          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);
+    }
   }
 }

Common/src/geometry/dual_grid/CPoint.cpp

@@ -132,6 +132,9 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) {
 
   RoughnessHeight.resize(npoint) = su2double(0.0);
   SharpEdge_Distance.resize(npoint) = su2double(0.0);
+
+  Normals.resize(npoint, 3) = su2double(0.0);
+
 }
 
 void CPoint::SetElems(const vector<vector<long> >& elemsMatrix) { Elem = CCompressedSparsePatternL(elemsMatrix); }