[WIP]: Flamelet progress variable (FPV) model

.github/workflows/regression.yml

@@ -60,12 +60,14 @@ jobs:
     strategy:
       fail-fast: false
       matrix: 
-        testscript: ['tutorials.py', 'parallel_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py', 'hybrid_regression_AD.py']
+        testscript: ['tutorials.py', 'parallel_regression.py', 'flamelet_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py']
         include:
           - testscript: 'tutorials.py'
             tag: MPI
           - testscript: 'parallel_regression.py'
             tag: MPI
+          - testscript: 'flamelet_regression.py'
+            tag: MPI
           - testscript: 'parallel_regression_AD.py'
             tag: MPI
           - testscript: 'serial_regression.py'
@@ -90,7 +92,7 @@ jobs:
         uses: docker://su2code/test-su2:20200303
         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_flamelet -s ${{matrix.testscript}}
   unit_tests:
     runs-on: ubuntu-latest
     name: Unit Tests

.gitignore

@@ -95,5 +95,10 @@ build/
 # ninja binary (build system)
 ninja
 
+<<<<<<< HEAD
+# vscode config files
+.vscode/*
+=======
 # Ignore vscode folder
 .vscode/
+>>>>>>> develop

Common/include/grid_movement/CSurfaceMovement.hpp

@@ -242,6 +242,15 @@ class CSurfaceMovement : public CGridMovement {
    */
   void CheckFFDDimension(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox, unsigned short iFFDBox);
 
+  /*!
+   * \brief Check the geometrical constraints of the FFD
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] FFDBox - Array with all the free forms FFDBoxes of the computation.
+   * \param[in] iFFDBox - Index of FFD box.
+   */
+  void CheckGeomConstr(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox, unsigned short iFFDBox);
+
   /*!
    * \brief Set the Parametric coordinates.
    * \param[in] geometry - Geometrical definition of the problem.

Common/include/linear_algebra/CSysVector.hpp

@@ -359,6 +359,15 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     for (auto iVar = 0ul; iVar < nVar; iVar++) vec_val[iPoint * nVar + iVar] = 0.0;
   }
 
+  /*!
+   * \brief Set the value of one variable to zero for one block.
+   * \param[in] iPoint - Index of the block being set to zero.
+   * \param[in] iVar - Index of the variable being set to zero.
+   */
+  inline void SetBlock_Zero(unsigned long iPoint, unsigned long iVar) {
+    vec_val[iPoint * nVar + iVar] = 0.0;
+  }
+
   /*!
    * \brief Set "block" to the vector.
    * \note Template param Overwrite can be set to false to update existing values.

Common/include/option_structure.hpp

@@ -76,7 +76,7 @@ const unsigned int MAX_PARAMETERS = 10;       /*!< \brief Maximum number of para
 const unsigned int MAX_NUMBER_PERIODIC = 10;  /*!< \brief Maximum number of periodic boundary conditions. */
 const unsigned int MAX_STRING_SIZE = 200;     /*!< \brief Maximum number of domains. */
 const unsigned int MAX_NUMBER_FFD = 15;       /*!< \brief Maximum number of FFDBoxes for the FFD. */
-enum: unsigned int{MAX_SOLS = 12};            /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */
+enum: unsigned int{MAX_SOLS = 13};            /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */
 const unsigned int MAX_TERMS = 6;             /*!< \brief Maximum number of terms in the numerical equations (dimension of solver container array). */
 const unsigned int MAX_ZONES = 3;             /*!< \brief Maximum number of zones. */
 const unsigned int MAX_FE_KINDS = 4;          /*!< \brief Maximum number of Finite Elements. */
@@ -481,6 +481,8 @@ enum RUNTIME_TYPE {
   RUNTIME_TRANS_SYS = 22,     /*!< \brief One-physics case, the code is solving the turbulence model. */
   RUNTIME_RADIATION_SYS = 23, /*!< \brief One-physics case, the code is solving the radiation model. */
   RUNTIME_ADJRAD_SYS = 24,    /*!< \brief One-physics case, the code is solving the adjoint radiation model. */
+  RUNTIME_SCALAR_SYS = 25,    /*!< \brief One-physics case, the code is solving the scalar model. */
+  RUNTIME_ADJSCALAR_SYS = 26, /*!< \brief One-physics case, the code is solving the adjoint scalar model. */
 };
 
 const int FLOW_SOL = 0;     /*!< \brief Position of the mean flow solution in the solver container array. */
@@ -498,6 +500,9 @@ const int ADJRAD_SOL = 8;   /*!< \brief Position of the continuous adjoint turbu
 const int MESH_SOL = 9;      /*!< \brief Position of the mesh solver. */
 const int ADJMESH_SOL = 10;   /*!< \brief Position of the adjoint of the mesh solver. */
 
+const int SCALAR_SOL = 11;     /*!< \brief Position of the scalar transport solution in the solver container array. */
+const int ADJSCALAR_SOL = 12;  /*!< \brief Position of the adjoint scalar transport solution in the solver container array. */
+
 const int FEA_SOL = 0;      /*!< \brief Position of the FEA equation in the solution solver array. */
 const int ADJFEA_SOL = 1;   /*!< \brief Position of the FEA adjoint equation in the solution solver array. */
 
@@ -549,15 +554,16 @@ static const MapType<std::string, ENUM_SPACE> Space_Map = {
  * \brief Types of fluid model
  */
 enum ENUM_FLUIDMODEL {
-  STANDARD_AIR = 0,       /*!< \brief Standard air gas model. */
-  IDEAL_GAS = 1,          /*!< \brief Ideal gas model. */
-  VW_GAS = 2,             /*!< \brief Van Der Waals gas model. */
-  PR_GAS = 3,             /*!< \brief Perfect Real gas model. */
-  CONSTANT_DENSITY = 4,   /*!< \brief Constant density gas model. */
-  INC_IDEAL_GAS = 5,      /*!< \brief Incompressible ideal gas model. */
-  INC_IDEAL_GAS_POLY = 6, /*!< \brief Inc. ideal gas, polynomial gas model. */
-  MUTATIONPP = 7,         /*!< \brief Mutation++ gas model for nonequilibrium flow. */
-  SU2_NONEQ = 8           /*!< \brief User defined gas model for nonequilibrium flow. */
+  STANDARD_AIR         = 0, /*!< \brief Standard air gas model. */
+  IDEAL_GAS            = 1, /*!< \brief Ideal gas model. */
+  VW_GAS               = 2, /*!< \brief Van Der Waals gas model. */
+  PR_GAS               = 3, /*!< \brief Perfect Real gas model. */
+  CONSTANT_DENSITY     = 4, /*!< \brief Constant density gas model. */
+  INC_IDEAL_GAS        = 5, /*!< \brief Incompressible ideal gas model. */
+  INC_IDEAL_GAS_POLY   = 6, /*!< \brief Inc. ideal gas, polynomial gas model. */
+  MUTATIONPP           = 7, /*!< \brief Mutation++ gas model for nonequilibrium flow. */
+  SU2_NONEQ            = 8, /*!< \brief User defined gas model for nonequilibrium flow. */
+  FLAMELET_FLUID_MODEL = 9, /*!, \brief Flamelet model */
 };
 static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("STANDARD_AIR", STANDARD_AIR)
@@ -568,6 +574,7 @@ static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("INC_IDEAL_GAS", INC_IDEAL_GAS)
   MakePair("INC_IDEAL_GAS_POLY", INC_IDEAL_GAS_POLY)
   MakePair("MUTATIONPP", MUTATIONPP)
+  MakePair("FLAMELET_FLUID_MODEL", FLAMELET_FLUID_MODEL)
   MakePair("SU2_NONEQ", SU2_NONEQ)
 };
 
@@ -652,34 +659,43 @@ static const MapType<std::string, FREESTREAM_OPTION> FreeStreamOption_Map = {
 /*!
  * \brief Types of viscosity model
  */
+
 enum class VISCOSITYMODEL {
   CONSTANT, /*!< \brief Constant viscosity. */
   SUTHERLAND, /*!< \brief Sutherlands Law viscosity. */
   POLYNOMIAL, /*!< \brief Polynomial viscosity. */
+  FLAMELET, /*!< \brief flamelet combustion model */
 };
 static const MapType<std::string, VISCOSITYMODEL> ViscosityModel_Map = {
   MakePair("CONSTANT_VISCOSITY", VISCOSITYMODEL::CONSTANT)
   MakePair("SUTHERLAND", VISCOSITYMODEL::SUTHERLAND)
   MakePair("POLYNOMIAL_VISCOSITY", VISCOSITYMODEL::POLYNOMIAL)
+  MakePair("FLAMELET", VISCOSITYMODEL::FLAMELET)
+
 };
 
 /*!
  * \brief Types of thermal conductivity model
  */
+
 enum class CONDUCTIVITYMODEL {
   CONSTANT, /*!< \brief Constant thermal conductivity. */
   CONSTANT_PRANDTL, /*!< \brief Constant Prandtl number. */
   POLYNOMIAL, /*!< \brief Polynomial thermal conductivity. */
+  FLAMELET, /*!< \brief flamelet combustion model */ 
 };
 static const MapType<std::string, CONDUCTIVITYMODEL> ConductivityModel_Map = {
   MakePair("CONSTANT_CONDUCTIVITY", CONDUCTIVITYMODEL::CONSTANT)
   MakePair("CONSTANT_PRANDTL", CONDUCTIVITYMODEL::CONSTANT_PRANDTL)
   MakePair("POLYNOMIAL_CONDUCTIVITY", CONDUCTIVITYMODEL::POLYNOMIAL)
+  MakePair("FLAMELET", CONDUCTIVITYMODEL::FLAMELET)
+
 };
 
 /*!
  * \brief Types of turbulent thermal conductivity model
  */
+
 enum class CONDUCTIVITYMODEL_TURB {
   NONE, /*!< \brief No turbulent contribution to the effective thermal conductivity for RANS. */
   CONSTANT_PRANDTL, /*!< \brief Include contribution to effective conductivity using constant turbulent Prandtl number for RANS. */
@@ -690,7 +706,22 @@ static const MapType<std::string, CONDUCTIVITYMODEL_TURB> TurbConductivityModel_
 };
 
 /*!
- * \brief Types of unsteady mesh motion
+ * \brief types of mass diffusivity models
+ */
+enum class DIFFUSIVITYMODEL {
+  CONSTANT_DIFFUSIVITY, /*!< \brief Constant mass diffusivity for scalar transport. */
+  CONSTANT_SCHMIDT,     /*!< \brief Constant Schmidt number for mass diffusion in scalar transport. */
+  FLAMELET,             /*!< \brief flamelet model */
+};
+
+static const MapType<std::string, DIFFUSIVITYMODEL> DiffusivityModel_Map = {
+  MakePair("CONSTANT_DIFFUSIVITY", DIFFUSIVITYMODEL::CONSTANT_DIFFUSIVITY)
+  MakePair("CONSTANT_SCHMIDT", DIFFUSIVITYMODEL::CONSTANT_SCHMIDT)
+  MakePair("FLAMELET", DIFFUSIVITYMODEL::FLAMELET)
+};
+
+/*!
+ * \brief types of unsteady mesh motion
  */
 enum ENUM_GRIDMOVEMENT {
   NO_MOVEMENT = 0,          /*!< \brief Simulation on a static mesh. */
@@ -737,7 +768,7 @@ enum ENUM_GUST_TYPE {
   SINE = 2,         /*!< \brief Sine shaped gust */
   ONE_M_COSINE = 3, /*!< \brief 1-cosine shaped gust */
   VORTEX = 4,       /*!< \brief A gust made from vortices */
-  EOG = 5           /*!< \brief An extreme operating gust */
+  EOG = 5,          /*!< \brief An extreme operating gust */
 };
 static const MapType<std::string, ENUM_GUST_TYPE> Gust_Type_Map = {
   MakePair("NONE", NO_GUST)
@@ -924,7 +955,23 @@ static const MapType<std::string, ENUM_TRANS_MODEL> Trans_Model_Map = {
 };
 
 /*!
- * \brief Types of subgrid scale models
+ * \brief types of scalar transport models
+ */
+enum ENUM_SCALAR_MODEL {
+  NO_SCALAR_MODEL   = 0,  /*!< \brief No scalar transport model. */
+  PASSIVE_SCALAR    = 1,  /*!< \brief Passive scalar transport model. */
+  PROGRESS_VARIABLE = 2,  /*!< \brief Progress variable combustion model. */
+  CUSTOM_SCALAR     = 3,  /*!< \brief Custom scalar transport model. */
+};
+static const MapType<std::string, ENUM_SCALAR_MODEL> Scalar_Model_Map = {
+  MakePair("NONE", NO_SCALAR_MODEL)
+  MakePair("PASSIVE_SCALAR", PASSIVE_SCALAR)
+  MakePair("PROGRESS_VARIABLE", PROGRESS_VARIABLE)
+  MakePair("CUSTOM_SCALAR", CUSTOM_SCALAR)
+};
+
+/*!
+ * \brief types of subgrid scale models
  */
 enum ENUM_SGS_MODEL {
   NO_SGS_MODEL = 0, /*!< \brief No subgrid scale model. */
@@ -1523,7 +1570,11 @@ enum ENUM_OBJECTIVE {
   TOPOL_DISCRETENESS = 63,      /*!< \brief Measure of the discreteness of the current topology. */
   TOPOL_COMPLIANCE = 64,        /*!< \brief Measure of the discreteness of the current topology. */
   STRESS_PENALTY = 65,          /*!< \brief Penalty function of VM stresses above a maximum value. */
+  SURFACE_CO = 80,
+  SURFACE_NOX = 81,
+  SURFACE_TEMP = 82,
 };
+
 static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("DRAG", DRAG_COEFFICIENT)
   MakePair("LIFT", LIFT_COEFFICIENT)
@@ -1575,6 +1626,9 @@ static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("TOPOL_DISCRETENESS", TOPOL_DISCRETENESS)
   MakePair("TOPOL_COMPLIANCE", TOPOL_COMPLIANCE)
   MakePair("STRESS_PENALTY", STRESS_PENALTY)
+  MakePair("AVG_CO", SURFACE_CO)
+  MakePair("AVG_NOX", SURFACE_NOX)
+  MakePair("AVG_TEMP", SURFACE_TEMP)
 };
 
 /*!
@@ -2249,6 +2303,20 @@ static const MapType<std::string, ENUM_VERIFICATION_SOLUTIONS> Verification_Solu
   MakePair("USER_DEFINED_SOLUTION",    USER_DEFINED_SOLUTION)
 };
 
+  /* the order matters: */
+  /* enthalpy and progress variable are the first and second entry in the lookup table */
+  enum FLAMELET_SCALAR_VARIABLES {
+    I_PROG_VAR,
+    I_ENTHALPY,
+    I_CO,
+    I_NOX,
+  };
+
+  enum FLAMELET_SCALAR_SOURCES {
+    I_SRC_TOT_PROG_VAR,
+    I_SRC_TOT_CO,
+    I_SRC_TOT_NOX,
+  };
 /*!
  * \brief Types of streamwise periodicity.
  */

Common/include/option_structure.inl

@@ -1101,6 +1101,78 @@ public:
   }
 };
 
+class COptionInletScalar : public COptionBase {
+  string name; // identifier for the option
+  unsigned short & size;
+  string * & marker;
+  su2double ** & inletscalarval; 
+
+public:
+  COptionInletScalar(string option_field_name, unsigned short & nMarker_Inlet_Scalar, string* & Marker_Inlet_Scalar, su2double** & option_field) : size(nMarker_Inlet_Scalar), marker(Marker_Inlet_Scalar), inletscalarval(option_field) {
+    this->name = option_field_name;
+  }
+
+  ~COptionInletScalar() {};
+
+  string SetValue(vector<string> option_value) {
+    COptionBase::SetValue(option_value);
+    unsigned short option_size = option_value.size();
+    if ((option_size == 1) && (option_value[0].compare("NONE") == 0)) {
+      this->size = 0;
+      this->marker = NULL;
+      this->inletscalarval = NULL;     
+      return "";
+    }
+
+    // at this point we need to know how many inlets we have. We count the number of strings
+    // in the list. We assume that the name of the inlet does not start with a number
+    unsigned short nVals=0;
+    unsigned short nScalars=0;
+    for (unsigned long i = 0; i < option_size; i++) {
+      // get the next entry in the inlet list
+      if (isalpha(option_value[i][0])) {
+        nVals++;
+      } else {
+        nScalars++;
+      }
+
+    }  
+
+
+    nScalars = nScalars/nVals;
+
+
+    // nVals is the number of inlets
+    this->size = nVals;
+    this->marker = new string[nVals];
+    this->inletscalarval = new su2double*[nVals];
+    for (unsigned long i = 0; i < nVals; i++) {
+
+      this->inletscalarval[i] = new su2double[nScalars];
+    }    
+
+    for (unsigned long i = 0; i < nVals; i++) {
+      this->marker[i].assign(option_value[(1+nScalars)*i]);
+
+      for (unsigned long j=0;j<nScalars;j++){
+        istringstream ss_nd(option_value[(1+nScalars)*i + 1+j]);
+        if (!(ss_nd >> this->inletscalarval[i][j])) {
+          return badValue(option_value, "inlet", this->name);
+        }
+      }
+
+    }
+
+    return "";
+  }
+
+  void SetDefault() {
+    this->marker = NULL;    
+    this->inletscalarval = NULL;
+    this->size = 0; // There is no default value for list
+  }
+};
+
 template <class Tenum>
 class COptionRiemann : public COptionBase {