[WIP] Implementation of Amplification Factor Transport(AFT) 2019b tranistion model

Common/include/CConfig.hpp

@@ -628,6 +628,7 @@ class CConfig {
   TURB_SGS_MODEL Kind_SGS_Model;    /*!< \brief LES SGS model definition. */
   TURB_TRANS_MODEL Kind_Trans_Model;  /*!< \brief Transition model definition. */
   TURB_TRANS_CORRELATION Kind_Trans_Correlation;  /*!< \brief Transition correlation model definition. */
+  AFT_CORRELATION Kind_AFT_Correlation;  /*!< \brief Transition correlation of AFT model definition. */
   su2double hRoughness;             /*!< \brief RMS roughness for Transition model. */
   unsigned short Kind_ActDisk, Kind_Engine_Inflow,
   *Kind_Data_Riemann,
@@ -765,11 +766,13 @@ class CConfig {
   string *Config_Filenames;           /*!< \brief List of names for configuration files. */
   SST_OPTIONS *SST_Options;           /*!< \brief List of modifications/corrections/versions of SST turbulence model.*/
   SA_OPTIONS *SA_Options;             /*!< \brief List of modifications/corrections/versions of SA turbulence model.*/
-  LM_OPTIONS *LM_Options;             /*!< \brief List of modifications/corrections/versions of SA turbulence model.*/
+  LM_OPTIONS *LM_Options;             /*!< \brief List of modifications/corrections/versions of LM transition model.*/
+  AFT_OPTIONS *AFT_Options;           /*!< \brief List of modifications/corrections/versions of AFT transition model.*/
   ROUGHSST_MODEL Kind_RoughSST_Model;         /*!< \brief List of modifications/corrections/versions of rough-wall boundary conditions for SST turbulence model.*/
   unsigned short nSST_Options;        /*!< \brief Number of SST options specified. */
   unsigned short nSA_Options;         /*!< \brief Number of SA options specified. */
-  unsigned short nLM_Options;         /*!< \brief Number of SA options specified. */
+  unsigned short nLM_Options;         /*!< \brief Number of LM options specified. */
+  unsigned short nAFT_Options;        /*!< \brief Number of AFT options specified. */
   WALL_FUNCTIONS  *Kind_WallFunctions;        /*!< \brief The kind of wall function to use for the corresponding markers. */
   unsigned short  **IntInfo_WallFunctions;    /*!< \brief Additional integer information for the wall function markers. */
   su2double       **DoubleInfo_WallFunctions; /*!< \brief Additional double information for the wall function markers. */
@@ -922,6 +925,7 @@ class CConfig {
   Tke_FreeStream,                  /*!< \brief Total turbulent kinetic energy of the fluid.  */
   Intermittency_FreeStream,        /*!< \brief Freestream intermittency (for sagt transition model) of the fluid.  */
   ReThetaT_FreeStream,             /*!< \brief Freestream Transition Momentum Thickness Reynolds Number (for LM transition model) of the fluid.  */
+  N_Critcal,                       /*!< \brief Critical N-factor (for AFT model).  */
   NuFactor_FreeStream,             /*!< \brief Ratio of turbulent to laminar viscosity. */
   NuFactor_Engine,                 /*!< \brief Ratio of turbulent to laminar viscosity at the engine. */
   KFactor_LowerLimit,               /*!< \Non dimensional coefficient for lower limit of K in SST model. */
@@ -1237,6 +1241,7 @@ class CConfig {
   SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */
   SA_ParsedOptions saParsedOptions;   /*!< \brief Additional parameters for the SA turbulence model. */
   LM_ParsedOptions lmParsedOptions;   /*!< \brief Additional parameters for the LM transition model. */
+  AFT_ParsedOptions aftParsedOptions; /*!< \brief Additional parameters for the AFT transition model. */
   su2double uq_delta_b;         /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */
   unsigned short eig_val_comp;  /*!< \brief Parameter used to determine type of eigenvalue perturbation */
   su2double uq_urlx;            /*!< \brief Under-relaxation factor */
@@ -2088,6 +2093,12 @@ class CConfig {
    */
   su2double GetReThetaT_FreeStream() const { return ReThetaT_FreeStream; }
 
+  /*!
+   * \brief Get the value of the critical N-factor.
+   * \return The critical N-factor.
+   */
+  su2double GetN_Critical(void) const { return N_Critcal; }
+
   /*!
    * \brief Get the value of the non-dimensionalized freestream turbulence intensity.
    * \return Non-dimensionalized freestream intensity.
@@ -2794,6 +2805,12 @@ class CConfig {
    */
   void SetReThetaT_FreeStream(su2double val_ReThetaT_freestream) { ReThetaT_FreeStream = val_ReThetaT_freestream; }
 
+  /*!
+   * \brief Set the freestream momentum thickness Reynolds number.
+   * \param[in] val_ReThetaT_freestream - Value of the freestream momentum thickness Reynolds number.
+   */
+  void SetN_Crtical(su2double val_N_critcal) { N_Critcal = val_N_critcal; }
+
   /*!
    * \brief Set the non-dimensional freestream energy.
    * \param[in] val_energy_freestreamnd - Value of the non-dimensional freestream energy.
@@ -10215,6 +10232,12 @@ class CConfig {
    */
   LM_ParsedOptions GetLMParsedOptions() const { return lmParsedOptions; }
 
+  /*!
+   * \brief Get parsed AFT option data structure.
+   * \return AFT option data structure.
+   */
+  AFT_ParsedOptions GetAFTParsedOptions() const { return aftParsedOptions; }
+
   /*!
    * \brief Get rough-wall boundary conditions for SST.
    */

Common/include/adt/CADTElemClass.hpp

@@ -122,6 +122,25 @@ class CADTElemClass : public CADTBaseClass {
                                  markerID, elemID, rankID);
   }
 
+  /*!
+   * \brief Determines the nearest element in the ADT for a given coordinate and computes
+   *        the closest point on that element.
+   * \note This function forwards the call to the implementation routine, selecting the
+   *       appropriate thread-local working variables.
+   * \param[in]  coor     Coordinate for which the nearest element must be determined.
+   * \param[out] dist     Distance to the nearest element.
+   * \param[out] xw       Coordinates of the closest point on the nearest element.
+   * \param[out] markerID Local marker ID of the nearest element.
+   * \param[out] elemID   Local element ID of the nearest element.
+   * \param[out] rankID   Rank on which the nearest element is stored.
+   */ 
+  inline void DetermineNearestElement(const su2double* coor, su2double& dist, su2double* xw,
+                                            unsigned short& markerID, unsigned long& elemID, int& rankID) {
+    const auto iThread = omp_get_thread_num();
+    DetermineNearestElement_impl(BBoxTargets[iThread], FrontLeaves[iThread], FrontLeavesNew[iThread], coor, dist,
+                                 xw, markerID, elemID, rankID);
+  }
+
  private:
   /*!
    * \brief Implementation of DetermineContainingElement.
@@ -139,6 +158,17 @@ class CADTElemClass : public CADTBaseClass {
                                     vector<unsigned long>& frontLeavesNew, const su2double* coor, su2double& dist,
                                     unsigned short& markerID, unsigned long& elemID, int& rankID) const;
 
+  /*!
+   * \brief Implementation of the nearest element search in the ADT for a given coordinate,
+   *        including computation of the closest point on the element.
+   * \note This routine performs the actual traversal of the ADT using the provided
+   *       thread-local working arrays. It identifies the nearest element and evaluates
+   *       the closest point on that element surface.
+   */
+  void DetermineNearestElement_impl(vector<CBBoxTargetClass>& BBoxTargets, vector<unsigned long>& frontLeaves,
+                                    vector<unsigned long>& frontLeavesNew, const su2double* coor, su2double& dist,
+                                    su2double* xw, unsigned short& markerID, unsigned long& elemID, int& rankID) const;
+
   /*!
    * \brief Function, which checks whether or not the given coordinate is
             inside the given element.
@@ -292,6 +322,15 @@ class CADTElemClass : public CADTBaseClass {
    * \param[out] dist2Elem Distance squared from the coordinate to the element.
    */
   void Dist2ToElement(const unsigned long elemID, const su2double* coor, su2double& dist2Elem) const;
+  /*!
+   * \brief Function, which computes the distance squared of the given coordinate
+            to the given element.
+   * \param[in]  elemID    ID of the element to which the distance must be determined.
+   * \param[in]  coor      Coordinate for which the distance to the element must be determined.
+   * \param[out] dist2Elem Distance squared from the coordinate to the element.
+   * \param[out] xw        Representative point used for the final distance evaluation.
+   */
+  void Dist2ToElement(const unsigned long elemID, const su2double* coor, su2double& dist2Elem, su2double* xw) const;
   /*!
    * \brief Function, which computes the distance squared of the given coordinate
             to a linear line element.
@@ -303,6 +342,19 @@ class CADTElemClass : public CADTBaseClass {
    * \param[out] dist2Line Distance squared from the coordinate to the line.
    */
   void Dist2ToLine(const unsigned long i0, const unsigned long i1, const su2double* coor, su2double& dist2Line) const;
+  /*!
+   * \brief Function, which computes the distance squared of the given coordinate to a linear
+   *        line element and returns the representative point used for the final distance evaluation.
+   * \param[in]  i0        Starting index in coorPoints, where the coordinates of the first point
+   *                       of the line are stored.
+   * \param[in]  i1        Starting index in coorPoints, where the coordinates of the second point
+   *                       of the line are stored.
+   * \param[in]  coor      Coordinate for which the distance to the line must be determined.
+   * \param[out] dist2Line Distance squared from the coordinate to the line.
+   * \param[out] xw        Representative point used for the final distance evaluation.
+   */
+  void Dist2ToLine(const unsigned long i0, const unsigned long i1, const su2double* coor,
+                 su2double& dist2Line, su2double* xw) const;
   /*!
    * \brief Function, which computes the distance squared of the given coordinate
             to a linear quadrilateral element if the projection is inside the quad.

Common/include/geometry/CGeometry.hpp

@@ -143,6 +143,8 @@ class CGeometry {
 
   vector<su2double> XCoordList; /*!< \brief Vector containing points appearing on a single plane */
 
+  su2activematrix WallNormalUnitVector; /*!< \brief Store wall-normal unit vector */
+
 #if defined(HAVE_MPI) && defined(HAVE_PARMETIS)
   vector<vector<unsigned long>>
       adj_nodes; /*!< \brief Vector of vectors holding each node's adjacency during preparation for ParMETIS. */
@@ -1887,4 +1889,29 @@ class CGeometry {
    * \return A pointer to the reference node coordinate vector.
    */
   inline virtual const su2double* GetStreamwise_Periodic_RefNode() const { return nullptr; }
+
+  public:
+  inline void AllocateWallNormalUnitVector() {
+    WallNormalUnitVector.resize(nPoint, nDim) = su2double(0.0);
+  }
+
+  /*!
+   * \brief Get a pointer to the reference node coordinate vector.
+   * \param iPoint - Point
+   * \param iDim - Dimension
+   * \param val - Wall-normal unit vector
+   */
+  inline void SetWallNormalUnitVector(unsigned long iPoint, unsigned short iDim, su2double val) {
+    WallNormalUnitVector(iPoint, iDim) = val;
+  }
+
+  inline su2double GetWallNormalUnitVector(unsigned long iPoint, unsigned short iDim) const {
+    return WallNormalUnitVector(iPoint, iDim);
+  }
+
+  inline const su2double* GetWallNormalUnitVector(unsigned long iPoint) const {
+  return &WallNormalUnitVector(iPoint, 0);
+  }
+
+  inline const auto& GetWallNormalUnitVector() const { return WallNormalUnitVector; }
 };

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. */
@@ -549,6 +550,26 @@ class CPoint {
    */
   inline su2double GetMaxLength(unsigned long iPoint) const { return MaxLength(iPoint); }
 
+  /*!
+   * \brief Set the value of the normal of the nearest wall element.
+   * \param[in] iPoint - Index of the point.
+   * \param[in] distance - Value of the normal.
+   */
+  template<typename Normals_type>
+  inline void SetNormal(unsigned long iPoint, Normals_type const&normalwall) { 
+    for (unsigned long iDim = 0; iDim <nDim; iDim++) {
+      Normals(iPoint, iDim) = normalwall[iDim];
+    } 
+  }
+
+  /*!
+   * \brief Get 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 Get area or volume of the control volume.
    * \param[in] iPoint - Index of the point.
@@ -885,4 +906,20 @@ class CPoint {
       }
     }
   }
-};
+
+  /*!
+   * \brief Set wall normal according to stored closest wall information.
+   * \param[in] normalwall - Mapping [rank][zone][marker][element] -> normalwall
+   */
+  template<typename Normals_type>
+  void SetWallNormals(Normals_type const&normalwall){
+    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, normalwall[rankID][zoneID][markerID][elementID]);
+    }
+  }
+};

Common/include/option_structure.hpp

@@ -1222,10 +1222,12 @@ inline SA_ParsedOptions ParseSAOptions(const SA_OPTIONS *SA_Options, unsigned sh
 enum class TURB_TRANS_MODEL {
   NONE,  /*!< \brief No transition model. */
   LM,    /*!< \brief Kind of transition model (Langtry-Menter (LM) for SST and Spalart-Allmaras). */
+  AFT,   /*!< \brief Kind of transition model (Amplification Factor Transport model for Spalart-Allmaras). */
 };
 static const MapType<std::string, TURB_TRANS_MODEL> Trans_Model_Map = {
   MakePair("NONE", TURB_TRANS_MODEL::NONE)
   MakePair("LM", TURB_TRANS_MODEL::LM)
+  MakePair("AFT", TURB_TRANS_MODEL::AFT)
 };
 
 /*!
@@ -1283,10 +1285,10 @@ struct LM_ParsedOptions {
  * \brief Function to parse LM options.
  * \param[in] LM_Options - Selected LM option from config.
  * \param[in] nLM_Options - Number of options selected.
- * \param[in] rank - MPI rank.
- * \return Struct with SA options.
+ * \param[in] Kind_Turb_Model- Selected based turbulence model from config.
+ * \return Struct with LM options.
  */
-inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned short nLM_Options, int rank, TURB_MODEL Kind_Turb_Model) {
+inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned short nLM_Options, TURB_MODEL Kind_Turb_Model) {
   LM_ParsedOptions LMParsedOptions;
 
   auto IsPresent = [&](LM_OPTIONS option) {
@@ -1342,7 +1344,58 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh
 }
 
 /*!
- * \brief Structure containing parsed options for data-driven fluid model.
+ * \brief AFT Options
+ */
+enum class AFT_OPTIONS {
+  NONE,         /*!< \brief No option / default. */
+  AFT2019b      /*!< \brief 2019b Coder SA-AFT model (https://doi.org/10.2514/6.2019-0039). */
+};
+
+static const MapType<std::string, AFT_OPTIONS> AFT_Options_Map = {
+  MakePair("NONE", AFT_OPTIONS::NONE)
+  MakePair("AFT2019b", AFT_OPTIONS::AFT2019b)
+};
+
+/*!
+ * \brief Types of transition correlations
+ */
+enum class AFT_CORRELATION {
+  NONE,         /*!< \brief No option / default. */
+  AFT2019b      /*!< \brief Kind of transition correlation model (AFT2019b). */
+};
+
+/*!
+ * \brief Structure containing parsed AFT options.
+ */
+struct AFT_ParsedOptions {
+  AFT_OPTIONS version = AFT_OPTIONS::NONE;  /*!< \brief AFT base model. */
+  AFT_CORRELATION Correlation = AFT_CORRELATION::NONE;
+};
+
+/*!
+ * \brief Function to parse AFT options.
+ * \param[in] AFT_Options - Selected AFT option from config.
+ * \param[in] nAFT_Options - Number of options selected.
+ * \return Struct with AFT options.
+ */
+inline AFT_ParsedOptions ParseAFTOptions(const AFT_OPTIONS *AFT_Options, unsigned short nAFT_Options) {
+  AFT_ParsedOptions AFTParsedOptions;
+
+  auto IsPresent = [&](AFT_OPTIONS option) {
+    const auto aft_options_end = AFT_Options + nAFT_Options;
+    return std::find(AFT_Options, aft_options_end, option) != aft_options_end;
+  };
+
+  int NFoundCorrelations = 0;
+  if (IsPresent(AFT_OPTIONS::AFT2019b)) {
+    AFTParsedOptions.Correlation = AFT_CORRELATION::AFT2019b;
+    AFTParsedOptions.version = AFT_OPTIONS::AFT2019b;
+    NFoundCorrelations++;
+  }
+  return AFTParsedOptions;
+}
+
+ /* \brief Structure containing parsed options for data-driven fluid model.
  */
 struct DataDrivenFluid_ParsedOptions {
   su2double Newton_relaxation = 1.0;  /*!< \brief Relaxation factor for Newton solvers in data-driven fluid models. */
@@ -2055,6 +2108,10 @@ enum class WALL_TYPE {
   SMOOTH,  /*!< \brief Smooth wall */
   ROUGH,   /*!< \brief Rough wall */
 };
+static const MapType<std::string, WALL_TYPE> WallType_Map = {
+  MakePair("SMOOTH", WALL_TYPE::SMOOTH)
+  MakePair("ROUGH", WALL_TYPE::ROUGH)
+};
 
 /*!
  * \brief Types of objective functions