Merge pull request #2057 from su2code/feature_new_flamelet

Feature new flamelet: partial premixed flamelet-generated manifold simulations are enabled through these features.

Author: EvertBunschoten

Common/include/CConfig.hpp

@@ -909,6 +909,7 @@ class CConfig {
   Tke_FreeStreamND,           /*!< \brief Farfield kinetic energy (external flow). */
   Omega_FreeStreamND,         /*!< \brief Specific dissipation (external flow). */
   Omega_FreeStream;           /*!< \brief Specific dissipation (external flow). */
+  bool Variable_Density;      /*!< \brief Variable density for incompressible flow. */
   unsigned short nElectric_Constant;    /*!< \brief Number of different electric constants. */
   su2double *Electric_Constant;         /*!< \brief Dielectric constant modulus. */
   su2double Knowles_B,                  /*!< \brief Knowles material model constant B. */
@@ -1229,9 +1230,10 @@ class CConfig {
   unsigned short n_user_sources = 0;  /*!< \brief Number of source terms for user defined (auxiliary) scalar transport equations. */
   unsigned short n_control_vars = 0;  /*!< \brief Number of controlling variables (independent variables) for the LUT. */
 
+  string* controlling_variable_names;
+  string* cv_source_names;
   vector<string> table_scalar_names;  /*!< \brief Names of transported scalar variables. */
-  string* table_lookup_names;         /*!< \brief Names of LUT variables. */
-  string file_name_lut;               /*!< \brief Filename of the LUT. */
+  string* lookup_names;         /*!< \brief Names of passive look-up variables. */
   string* user_scalar_names;          /*!< \brief Names of the user defined (auxiliary) transported scalars .*/
   string* user_source_names;          /*!< \brief Names of the source terms for the user defined transported scalars. */
 
@@ -2151,6 +2153,19 @@ class CConfig {
    */
   unsigned short GetNUserScalars() const { return n_user_scalars; }
 
+  /*!
+   * \brief Get the name of a specific controlling variable.
+   */
+  const string& GetControllingVariableName(unsigned short i_cv) const {
+    return controlling_variable_names[i_cv];
+  }
+
+  /*!
+   * \brief Get the name of the source term variable for a specific controlling variable.
+   */
+  const string& GetControllingVariableSourceName(unsigned short i_cv) const {
+    return cv_source_names[i_cv];
+  }
   /*!
    * \brief Get the name of the user scalar.
    */
@@ -2175,13 +2190,7 @@ class CConfig {
   /*!
    * \brief Get the name of the variable that we want to retrieve from the lookup table.
    */
-  const string& GetLUTLookupName(unsigned short i_lookup) const { return table_lookup_names[i_lookup]; }
-
-  /*!
-   * \brief Get the file name of the look up table.
-   * \return File name of the look up table.
-   */
-  const string& GetFileNameLUT() const { return file_name_lut; }
+  const string& GetLookupName(unsigned short i_lookup) const { return lookup_names[i_lookup]; }
 
   /*!
    * \brief Get the Young's modulus of elasticity.
@@ -3884,6 +3893,12 @@ class CConfig {
    */
   INC_DENSITYMODEL GetKind_DensityModel() const { return Kind_DensityModel; }
 
+  /*!
+   * \brief Selection of variable density option for incompressible flows.
+   * \return Flag for variable density for incompressible flows.
+   */
+  bool GetVariable_Density_Model() const { return Variable_Density; }
+
   /*!
    * \brief Flag for whether to solve the energy equation for incompressible flows.
    * \return Flag for energy equation

Common/include/containers/CLookUpTable.hpp

@@ -76,6 +76,8 @@ class CLookUpTable {
    */
   su2vector<su2activematrix> table_data;
 
+  double memory_footprint_data = 0; /*!< \brief Memory footprint of the loaded table data. */
+
   /*! \brief
    * Holds all connectivity data stored in the table for each level. First index
    * addresses the variable while second index addresses the point.
@@ -216,7 +218,7 @@ class CLookUpTable {
    */
   void Linear_Interpolation(const su2double val_CV3, const unsigned long lower_level, const unsigned long upper_level,
                             std::vector<su2double>& lower_values, std::vector<su2double>& upper_values,
-                            std::vector<su2double*>& var_vals) const;
+                            std::vector<su2double>& var_vals) const;
 
   /*!
    * \brief Find the point on the hull (boundary of the table) that is closest to the point P(val_CV1,val_CV2).
@@ -345,7 +347,7 @@ class CLookUpTable {
    * \param[in] val_CV3 - Value of controlling variable 3.
    * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not.
    */
-  unsigned long LookUp_XYZ(const std::vector<std::string>& val_names_var, std::vector<su2double*>& val_vars,
+  unsigned long LookUp_XYZ(const std::vector<std::string>& val_names_var, std::vector<su2double>& val_vars,
                            su2double val_CV1, su2double val_CV2, su2double val_CV3 = 0);
 
   /*!

Common/include/containers/CTrapezoidalMap.hpp

@@ -54,12 +54,14 @@ class CTrapezoidalMap {
    * same band. Used to sort the edges */
   su2vector<std::vector<std::pair<su2double, unsigned long> > > y_edge_at_band_mid;
 
+  double memory_footprint = 0;
+
  public:
   CTrapezoidalMap() = default;
 
   CTrapezoidalMap(const su2double* samples_x, const su2double* samples_y, const unsigned long size,
                   const std::vector<std::array<unsigned long, 2> >& edges,
-                  const su2vector<std::vector<unsigned long> >& edge_to_triangle);
+                  const su2vector<std::vector<unsigned long> >& edge_to_triangle, bool display = false);
 
   /*!
    * \brief return the index to the triangle that contains the coordinates (val_x,val_y)
@@ -97,4 +99,10 @@ class CTrapezoidalMap {
   inline bool IsInsideHullX(su2double val_x) {
     return (val_x >= unique_bands_x.front()) && (val_x <= unique_bands_x.back());
   }
+
+  /*!
+   * \brief get memory footprint of trapezoidal map.
+   * \return - memory footprint in mega bytes.
+   */
+  double GetMemoryFootprint() const { return memory_footprint; }
 };

Common/include/option_structure.hpp

@@ -628,11 +628,13 @@ enum class INC_DENSITYMODEL {
   CONSTANT,   /*!< \brief Constant density. */
   BOUSSINESQ, /*!< \brief Boussinesq density model. */
   VARIABLE,   /*!< \brief Variable density model. */
+  FLAMELET,   /*!< \brief Density according to flamelet manifold. */
 };
 static const MapType<std::string, INC_DENSITYMODEL> DensityModel_Map = {
   MakePair("CONSTANT", INC_DENSITYMODEL::CONSTANT)
   MakePair("BOUSSINESQ", INC_DENSITYMODEL::BOUSSINESQ)
   MakePair("VARIABLE", INC_DENSITYMODEL::VARIABLE)
+  MakePair("FLAMELET", INC_DENSITYMODEL::FLAMELET)
 };
 
 /*!
@@ -1317,6 +1319,7 @@ static const MapType<std::string, SPECIES_MODEL> Species_Model_Map = {
 enum FLAMELET_SCALAR_VARIABLES {
   I_PROGVAR,
   I_ENTH,
+  I_MIXFRAC,
 };
 
 /*!
@@ -1326,6 +1329,15 @@ enum FLAMELET_SCALAR_SOURCES {
   I_SRC_TOT_PROGVAR
 };
 
+/*!
+ * \brief Look-up operations for the flamelet scalar solver.
+ */
+enum FLAMELET_LOOKUP_OPS {
+  TD,       /*!< \brief Thermochemical properties (temperature, density, diffusivity, etc.). */
+  SOURCES,  /*!< \brief Scalar source terms (controlling variables, passive species).*/
+  LOOKUP,   /*!< \brief Passive look-up variables specified in config. */
+};
+
 /*!
  * \brief Types of subgrid scale models
  */

Common/src/CConfig.cpp

@@ -2101,11 +2101,14 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Determine if we need to allocate memory to store the multizone residual. \n DEFAULT: true (temporarily) */
   addBoolOption("MULTIZONE_RESIDUAL", Multizone_Residual, false);
 
-  /*!\brief File name of the flamelet look up table.*/
-  addStringOption("FILENAME_LUT", file_name_lut, string("LUT"));
+  /* !\brief CONTROLLING_VARIABLE_NAMES \n DESCRIPTION: Names of the variables used as inputs for the data regression method in flamelet or data-driven fluid models. */
+  addStringListOption("CONTROLLING_VARIABLE_NAMES", n_control_vars, controlling_variable_names);
+
+  /* !\brief CONTROLLING_VARIABLE_SOURCE_NAMES \n DESCRIPTION: Names of the variables in the flamelet manifold corresponding to the source terms of the controlling variables. */
+  addStringListOption("CONTROLLING_VARIABLE_SOURCE_NAMES", n_control_vars, cv_source_names);
 
   /* DESCRIPTION: Names of the passive lookup variables for flamelet LUT */
-  addStringListOption("LOOKUP_NAMES", n_lookups, table_lookup_names);
+  addStringListOption("LOOKUP_NAMES", n_lookups, lookup_names);
 
   /* DESCRIPTION: Names of the user transport equations solved in the flamelet problem. */
   addStringListOption("USER_SCALAR_NAMES", n_user_scalars, user_scalar_names);
@@ -3840,6 +3843,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   SetDefaultIfEmpty(Prandtl_Turb, nPrandtl_Turb, Prandtl_Turb_Default);
   SetDefaultIfEmpty(Constant_Lewis_Number, nConstant_Lewis_Number, Lewis_Number_Default);
 
+  Variable_Density = ((Kind_DensityModel == INC_DENSITYMODEL::VARIABLE) || (Kind_DensityModel == INC_DENSITYMODEL::FLAMELET));
+
   /*--- Check whether inputs for FLUID_MIXTURE are correctly specified. ---*/
 
     if (Kind_FluidModel == FLUID_MIXTURE) {
@@ -3935,8 +3940,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
                        CURRENT_FUNCTION);
       }
 
-      if (Kind_DensityModel != INC_DENSITYMODEL::VARIABLE) {
-        SU2_MPI::Error("The use of FLUID_FLAMELET requires the INC_DENSITY_MODEL option to be VARIABLE",
+      if (!Variable_Density) {
+        SU2_MPI::Error("The use of FLUID_FLAMELET requires the INC_DENSITY_MODEL option to be VARIABLE or FLAMELET",
                        CURRENT_FUNCTION);
       }
 
@@ -5480,8 +5485,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
   /*--- Define some variables for flamelet model. ---*/
   if (Kind_Species_Model == SPECIES_MODEL::FLAMELET) {
-    /*--- The controlling variables are progress variable and total enthalpy ---*/
-    n_control_vars = 2;
+    /*--- The controlling variables are progress variable, total enthalpy, and optionally mixture fraction ---*/
+    //n_control_vars = nSpecies - n_user_scalars;
+    if (n_control_vars != (nSpecies - n_user_scalars))
+      SU2_MPI::Error("Number of initial species incompatbile with number of controlling variables and user scalars.", CURRENT_FUNCTION);
     /*--- We can have additional user defined transported scalars ---*/
     n_scalars = n_control_vars + n_user_scalars;
   }

Common/src/containers/CLookUpTable.cpp

@@ -73,9 +73,13 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in,
   trap_map_x_y.resize(n_table_levels);
   su2double startTime = SU2_MPI::Wtime();
   unsigned short barwidth = 65;
+  bool display_map_info = (n_table_levels < 2);
+  double tmap_memory_footprint = 0;
   for (auto i_level = 0ul; i_level < n_table_levels; i_level++) {
-    trap_map_x_y[i_level] = CTrapezoidalMap(GetDataP(name_CV1, i_level), GetDataP(name_CV2, i_level),
-                                            table_data[i_level].cols(), edges[i_level], edge_to_triangle[i_level]);
+    trap_map_x_y[i_level] =
+        CTrapezoidalMap(GetDataP(name_CV1, i_level), GetDataP(name_CV2, i_level), table_data[i_level].cols(),
+                        edges[i_level], edge_to_triangle[i_level], display_map_info);
+    tmap_memory_footprint += trap_map_x_y[i_level].GetMemoryFootprint();
     /* Display a progress bar to monitor table generation process */
     if (rank == MASTER_NODE) {
       su2double progress = su2double(i_level) / n_table_levels;
@@ -100,6 +104,8 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in,
       default:
         break;
     }
+    cout << "Trapezoidal map memory footprint: " << tmap_memory_footprint << " MB\n";
+    cout << "Table data memory footprint: " << memory_footprint_data << " MB\n" << endl;
   }
 
   ComputeInterpCoeffs();
@@ -133,7 +139,9 @@ void CLookUpTable::LoadTableRaw(const string& var_file_name_lut) {
     table_data[i_level] = file_reader.GetTableData(i_level);
     triangles[i_level] = file_reader.GetTriangles(i_level);
     hull[i_level] = file_reader.GetHull(i_level);
+    memory_footprint_data += n_points[i_level] * sizeof(su2double);
   }
+  memory_footprint_data /= 1e6;
 
   n_variables = file_reader.GetNVariables();
   version_lut = file_reader.GetVersionLUT();
@@ -164,7 +172,7 @@ void CLookUpTable::FindTableLimits(const string& name_cv1, const string& name_cv
   }
 
   if (table_dim == 3) {
-    limits_table_z = minmax_element(z_values_levels.data(), z_values_levels.data() + z_values_levels.cols());
+    limits_table_z = minmax_element(z_values_levels.data(), z_values_levels.data() + z_values_levels.size());
   }
 }
 
@@ -451,7 +459,7 @@ unsigned long CLookUpTable::LookUp_XYZ(const std::string& val_name_var, su2doubl
     return 1;
   }
 }
-unsigned long CLookUpTable::LookUp_XYZ(const std::vector<std::string>& val_names_var, std::vector<su2double*>& val_vars,
+unsigned long CLookUpTable::LookUp_XYZ(const std::vector<std::string>& val_names_var, std::vector<su2double>& val_vars,
                                        su2double val_CV1, su2double val_CV2, su2double val_CV3) {
   /*--- Perform quasi-3D interpolation for a vector of variables with names val_names_var
         on a query point with coordinates val_CV1, val_CV2, and val_CV3 ---*/
@@ -509,7 +517,7 @@ void CLookUpTable::Linear_Interpolation(const su2double val_CV3, const unsigned
 
 void CLookUpTable::Linear_Interpolation(const su2double val_CV3, const unsigned long lower_level,
                                         const unsigned long upper_level, std::vector<su2double>& lower_values,
-                                        std::vector<su2double>& upper_values, std::vector<su2double*>& var_vals) const {
+                                        std::vector<su2double>& upper_values, std::vector<su2double>& var_vals) const {
   /* Perform linear interpolation along the z-direction of the table for multiple variables */
 
   /* Retrieve constant z-values of inclusion levels */
@@ -522,7 +530,7 @@ void CLookUpTable::Linear_Interpolation(const su2double val_CV3, const unsigned
 
   /* Perform linear interpolation */
   for (size_t iVar = 0; iVar < var_vals.size(); iVar++) {
-    *var_vals[iVar] = lower_values[iVar] * factor_lower + upper_values[iVar] * factor_upper;
+    var_vals[iVar] = lower_values[iVar] * factor_lower + upper_values[iVar] * factor_upper;
   }
 }
 

Common/src/containers/CTrapezoidalMap.cpp

@@ -40,7 +40,7 @@ using namespace std;
  */
 CTrapezoidalMap::CTrapezoidalMap(const su2double* samples_x, const su2double* samples_y, const unsigned long size,
                                  vector<std::array<unsigned long, 2> > const& edges,
-                                 su2vector<vector<unsigned long> > const& val_edge_to_triangle) {
+                                 su2vector<vector<unsigned long> > const& val_edge_to_triangle, bool display) {
   int rank = SU2_MPI::GetRank();
   su2double startTime = SU2_MPI::Wtime();
 
@@ -148,11 +148,11 @@ CTrapezoidalMap::CTrapezoidalMap(const su2double* samples_x, const su2double* sa
     for (unsigned long j_y = 0; j_y < y_edge_at_band_mid[i_y].size(); j_y++)
       size_y_edge_at_band_mid += sizeof(su2double) / 1e6 + sizeof(unsigned long) / 1e6;
 
-  double size_total =
+  memory_footprint =
       size_unique_bands + size_edge_limits_x + size_edge_limits_y + size_edge_to_triangle + size_y_edge_at_band_mid;
 
   /* print size of trapezoidal map components to screen */
-  if (rank == MASTER_NODE) {
+  if ((rank == MASTER_NODE) && display) {
     cout << setfill(' ');
     cout << "\n" << endl;
     cout << "+------------------------------------------------------------------+\n";
@@ -171,7 +171,7 @@ CTrapezoidalMap::CTrapezoidalMap(const su2double* samples_x, const su2double* sa
          << " |" << endl;
     cout << "| Size of y_edge_at_band_mid in memory: " << setw(22) << size_y_edge_at_band_mid << " MB "
          << " |" << endl;
-    cout << "| Total:                                " << setw(22) << size_total << " MB "
+    cout << "| Total:                                " << setw(22) << memory_footprint << " MB "
          << " |" << endl;
     cout << "+------------------------------------------------------------------+" << endl;
     cout << "\n" << endl;
@@ -198,9 +198,14 @@ unsigned long CTrapezoidalMap::GetTriangle(su2double val_x, su2double val_y) {
   sort(triangles_edge_up.begin(), triangles_edge_up.end());
 
   /* The intersection of the faces to which upper or lower belongs is the face that both belong to. */
-  vector<unsigned long> triangle(1);
+  vector<unsigned long> triangle;
   set_intersection(triangles_edge_up.begin(), triangles_edge_up.end(), triangles_edge_low.begin(),
-                   triangles_edge_low.end(), triangle.begin());
+                   triangles_edge_low.end(), std::back_inserter(triangle));
+
+  /*--- We failed to find an intersection, so take the lower triangle inside the band enclosing the point---*/
+  if (triangle.size() < 1) {
+    triangle.resize(1, triangles_edge_low[0]);
+  }
 
   return triangle[0];
 }

SU2_CFD/include/fluid/CDataDrivenFluid.hpp

@@ -210,7 +210,7 @@ class CDataDrivenFluid final : public CFluidModel {
    * \brief Get fluid model extrapolation instance.
    * \return Query point lies outside fluid model data range.
    */
-  unsigned long GetExtrapolation() override { return outside_dataset; }
+  unsigned long GetExtrapolation() const override { return outside_dataset; }
 
   /*!
    * \brief Get number of Newton solver iterations.