Feature new flamelet: addition of preferential diffusion

.github/workflows/regression.yml

@@ -171,7 +171,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:230813-0103
         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_new_flamelet -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:230813-0103
         with:

Common/include/CConfig.hpp

@@ -1220,8 +1220,13 @@ class CConfig {
   unsigned short nSpecies_Init;    /*!< \brief Number of entries of SPECIES_INIT */
 
   /*--- Additional flamelet solver options ---*/
-  su2double flame_init[8];       /*!< \brief Initial solution parameters for flamelet solver.*/
-
+  FLAMELET_INIT_TYPE flame_init_type = FLAMELET_INIT_TYPE::NONE; /*!< \brief Method for solution ignition for flamelet problems. */
+  su2double flame_init[8];       /*!< \brief Flame front initialization parameters. */
+  su2double spark_init[6];       /*!< \brief Spark ignition initialization parameters. */
+  su2double* spark_reaction_rates; /*!< \brief Source terms for flamelet spark ignition option. */
+  unsigned short nspark;           /*!< \brief Number of source terms for spark initialization. */
+  bool preferential_diffusion = false;  /*!< \brief Preferential diffusion physics for flamelet solver.*/
+
   /*--- lookup table ---*/
   unsigned short n_scalars = 0;       /*!< \brief Number of transported scalars for flamelet LUT approach. */
   unsigned short n_lookups = 0;       /*!< \brief Number of lookup variables, for visualization only. */
@@ -2130,13 +2135,44 @@ class CConfig {
 
   /*!
    * \brief Get the flame initialization.
-   *        (x1,x2,x3) = flame offset.
-   *        (x4,x5,x6) = flame normal, separating unburnt from burnt.
+   *        (x1,x2,x3) = flame offset/spark center location.
+   *        (x4,x5,x6) = flame normal, separating unburnt from burnt or
+   *                     spark radius, spark start iteration, spark duration.
    *        (x7) = flame thickness, the length from unburnt to burnt conditions.
    *        (x8) = flame burnt thickness, the length to stay at burnt conditions.
-   * \return Flame initialization for the flamelet model.
+   * \return Ignition initialization parameters for the flamelet model.
+   */
+  const su2double* GetFlameInit() const { 
+    switch (flame_init_type)
+    {
+    case FLAMELET_INIT_TYPE::FLAME_FRONT:
+      return flame_init;
+      break;
+    case FLAMELET_INIT_TYPE::SPARK:
+      return spark_init;
+      break;
+    default:
+      return nullptr;
+      break;
+    }
+  }
+
+  /*!
+   * \brief Get species net reaction rates applied during spark ignition.
+   */
+  const su2double* GetSpark() const {
+    return spark_reaction_rates;
+  }
+
+  /*!
+   * \brief Preferential diffusion combustion problem.
+   */
+  bool GetPreferentialDiffusion() const { return preferential_diffusion; }
+
+  /*!
+   * \brief Define preferential diffusion combustion problem.
    */
-  const su2double* GetFlameInit() const { return flame_init; }
+  inline void SetPreferentialDiffusion(bool input) { preferential_diffusion = input; }
 
   /*!
    * \brief Get the number of control variables for flamelet model.
@@ -2182,6 +2218,7 @@ class CConfig {
     if (n_user_sources > 0) return user_source_names[i_user_source]; else return none;
   }
 
+  FLAMELET_INIT_TYPE GetFlameletInitType() const { return flame_init_type; }
   /*!
    * \brief Get the number of transported scalars for combustion.
    */

Common/include/containers/CLookUpTable.hpp

@@ -373,4 +373,9 @@ class CLookUpTable {
   inline std::pair<su2double*, su2double*> GetTableLimitsX(unsigned long i_level = 0) const {
     return limits_table_x[i_level];
   }
+
+  /*!
+   * \brief Check whether provided set of variables are included in the table.
+   */
+  bool CheckForVariables(const std::vector<std::string>& vars_to_check);
 };

Common/include/option_structure.hpp

@@ -1334,10 +1334,37 @@ enum FLAMELET_SCALAR_SOURCES {
  */
 enum FLAMELET_LOOKUP_OPS {
   TD,       /*!< \brief Thermochemical properties (temperature, density, diffusivity, etc.). */
+  PD,       /*!< \brief Preferential diffusion scalars. */
   SOURCES,  /*!< \brief Scalar source terms (controlling variables, passive species).*/
   LOOKUP,   /*!< \brief Passive look-up variables specified in config. */
 };
 
+/*!
+ * \brief the preferential diffusion scalar indices for the preferential diffusion model.
+ */
+enum FLAMELET_PREF_DIFF_SCALARS {
+  I_BETA_PROGVAR,
+  I_BETA_ENTH_THERMAL,
+  I_BETA_ENTH,
+  I_BETA_MIXFRAC,
+  N_BETA_TERMS,
+};
+
+/*!
+ * \brief Flame initialization options for the flamelet solver.
+ */
+enum class FLAMELET_INIT_TYPE {
+  FLAME_FRONT,  /*!< \brief Straight flame front. */
+  SPARK,        /*!< \brief Species reaction rate in a set location. */
+  NONE,         /*!< \brief No ignition, cold flow only. */
+};
+
+static const MapType<std::string, FLAMELET_INIT_TYPE> Flamelet_Init_Map = {
+  MakePair("NONE", FLAMELET_INIT_TYPE::NONE)
+  MakePair("FLAME_FRONT", FLAMELET_INIT_TYPE::FLAME_FRONT)
+  MakePair("SPARK", FLAMELET_INIT_TYPE::SPARK)
+};
+
 /*!
  * \brief Types of subgrid scale models
  */

Common/src/CConfig.cpp

@@ -972,7 +972,7 @@
   Species_Init           = nullptr;
   Species_Clipping_Min   = nullptr;
   Species_Clipping_Max   = nullptr;
-
+  spark_reaction_rates   = nullptr;
   /*--- Moving mesh pointers ---*/
 
   nKind_SurfaceMovement = 0;
@@ -1361,7 +1361,9 @@
   /*!\brief SPECIES_CLIPPING_MIN \n DESCRIPTION: Minimum values for scalar clipping \ingroup Config*/
   addDoubleListOption("SPECIES_CLIPPING_MIN", nSpecies_Clipping_Min, Species_Clipping_Min);
 
-  /*!\brief FLAME_INIT \n DESCRIPTION: flame initialization using the flamelet model \ingroup Config*/
+  /*!\brief FLAME_INIT_METHOD \n DESCRIPTION: Ignition method for flamelet solver \n DEFAULT: no ignition; cold flow only. */
+  addEnumOption("FLAME_INIT_METHOD", flame_init_type, Flamelet_Init_Map, FLAMELET_INIT_TYPE::NONE);
+  /*!\brief FLAME_INIT \n DESCRIPTION: flame front initialization using the flamelet model \ingroup Config*/
   /*--- flame offset (x,y,z) ---*/
   flame_init[0] = 0.0; flame_init[1] = 0.0; flame_init[2] = 0.0;
   /*--- flame normal (nx, ny, nz) ---*/
@@ -1370,6 +1372,13 @@
   flame_init[6] = 0.5e-3; flame_init[7] = 1.0;
   addDoubleArrayOption("FLAME_INIT", 8,flame_init);
 
+  /*!\brief SPARK_INIT \n DESCRIPTION: spark initialization using the flamelet model \ingroup Config*/
+  for (auto iSpark=0u; iSpark<6; ++iSpark) spark_init[iSpark]=0;
+  addDoubleArrayOption("SPARK_INIT", 6, spark_init);
+
+  /*!\brief SPARK_REACTION_RATES \n DESCRIPTION: Net source term values applied to species within spark area during spark ignition. \ingroup Config*/
+  addDoubleListOption("SPARK_REACTION_RATES", nspark, spark_reaction_rates);
+
   /*--- Options related to mass diffusivity and thereby the species solver. ---*/
 
   /*!\brief DIFFUSIVITY_MODEL\n DESCRIPTION: mass diffusivity model \n DEFAULT constant disffusivity \ingroup Config*/
@@ -9168,7 +9177,7 @@
     if (Marker_Isothermal[iMarker_Isothermal] == val_marker)
       return Isothermal_Temperature[iMarker_Isothermal];
 
-  return Isothermal_Temperature[0];
+  return Temperature_FreeStream;// Isothermal_Temperature[0];
 }
 
 su2double CConfig::GetWall_HeatFlux(const string& val_marker) const {

Common/src/containers/CLookUpTable.cpp

@@ -800,3 +800,15 @@ bool CLookUpTable::IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned l
 
   return (abs(area_tri - (area_0 + area_1 + area_2)) < area_tri * 1e-10);
 }
+
+bool CLookUpTable::CheckForVariables(const std::vector<std::string>& vars_to_check) {
+  bool variables_are_present{true};
+  for (auto iVar = 0u; iVar < vars_to_check.size(); iVar++) {
+    bool var_is_present{false};
+    for (auto jVar = 0u; jVar < names_var.size(); jVar++) {
+      if (vars_to_check[iVar].compare(names_var[jVar]) == 0) var_is_present = true;
+    }
+    if (!var_is_present) variables_are_present = false;
+  }
+  return variables_are_present;
+}

Docs/docmain.hpp

@@ -225,9 +225,3 @@
  * \brief Classes for explicit (done by the programmer) vectorization (SIMD) of computations.
  * \ingroup Toolboxes
  */
-
-/*!
- * \defgroup Multi-Layer Perceptrons (MLP)
- * \brief Data look up and interpolation via dense, feed-forward multi-layer perceptrons.
- * \ingroup Toolboxes
- */

SU2_CFD/include/fluid/CFluidFlamelet.hpp

@@ -42,9 +42,11 @@ class CFluidFlamelet final : public CFluidModel {
 
   enum LOOKUP_TD { TEMPERATURE, HEATCAPACITY, VISCOSITY, CONDUCTIVITY, DIFFUSIONCOEFFICIENT, MOLARWEIGHT, SIZE };
 
+  su2double beta_progvar, beta_enth_thermal, beta_enth, beta_mixfrac;
   int rank;
 
-  bool include_mixture_fraction = false; /*!< \brief include mixture fraction in controlling variables. */
+  bool include_mixture_fraction = false, /*!< \brief include mixture fraction in controlling variables. */
+      PreferentialDiffusion = false;     /*!< \brief use preferential diffusion physics. */
 
   unsigned short n_scalars, n_lookups, n_user_scalars, /*!< \brief number of passive reactant species. */
       n_control_vars;                                  /*!< \brief number of controlling variables. */
@@ -64,6 +66,7 @@ class CFluidFlamelet final : public CFluidModel {
 #ifdef USE_MLPCPP
   MLPToolbox::CLookUp_ANN* lookup_mlp; /*!< \brief Multi-layer perceptron collection. */
   MLPToolbox::CIOMap* iomap_TD;        /*!< \brief Input-output map for thermochemical properties. */
+  MLPToolbox::CIOMap* iomap_PD;        /*!< \brief Input-output map for the preferential diffusion scalars. */
   MLPToolbox::CIOMap* iomap_Sources;   /*!< \brief Input-output map for species source terms. */
   MLPToolbox::CIOMap* iomap_LookUp;    /*!< \brief Input-output map for passive look-up terms. */
   MLPToolbox::CIOMap* iomap_Current;
@@ -72,10 +75,10 @@ class CFluidFlamelet final : public CFluidModel {
   vector<su2double> scalars_vector;
 
   vector<string> varnames_TD, /*!< \brief Lookup names for thermodynamic state variables. */
-      varnames_Sources, varnames_LookUp;
+      varnames_Sources, varnames_LookUp, varnames_PD;
 
   vector<su2double> val_vars_TD, /*!< \brief References to thermodynamic state variables. */
-      val_vars_Sources, val_vars_LookUp;
+      val_vars_Sources, val_vars_LookUp, val_vars_PD;
 
   void PreprocessLookUp(CConfig* config);
 
@@ -127,4 +130,10 @@ class CFluidFlamelet final : public CFluidModel {
    * \param[out] Mu - value of the laminar viscosity
    */
   inline su2double GetLaminarViscosity() override { return Mu; }
+
+  /*!
+   * \brief Preferential diffusion as relevant phenomenon in flamelet simulations.
+   * \return Inclusion of preferential diffusion model.
+   */
+  inline bool GetPreferentialDiffusion() const override { return PreferentialDiffusion; };
 };

SU2_CFD/include/fluid/CFluidModel.hpp

@@ -46,34 +46,34 @@ class CLookUpTable;
  */
 class CFluidModel {
  protected:
-  su2double StaticEnergy{0.0}; /*!< \brief Internal Energy. */
-  su2double Entropy{0.0};      /*!< \brief Entropy. */
-  su2double Density{0.0};      /*!< \brief Density. */
-  su2double Pressure{0.0};     /*!< \brief Pressure. */
-  su2double SoundSpeed2{0.0};  /*!< \brief SpeedSound. */
-  su2double Temperature{0.0};  /*!< \brief Temperature. */
-  su2double dPdrho_e{0.0};     /*!< \brief DpDd_e. */
-  su2double dPde_rho{0.0};     /*!< \brief DpDe_d. */
-  su2double dTdrho_e{0.0};     /*!< \brief DTDd_e. */
-  su2double dTde_rho{0.0};     /*!< \brief DTDe_d. */
-  su2double dhdrho_P{0.0};     /*!< \brief DhDrho_p. */
-  su2double dhdP_rho{0.0};     /*!< \brief DhDp_rho. */
-  su2double dsdrho_P{0.0};     /*!< \brief DsDrho_p. */
-  su2double dsdP_rho{0.0};     /*!< \brief DsDp_rho. */
-  su2double Cp{0.0};           /*!< \brief Specific Heat Capacity at constant pressure. */
-  su2double Cv{0.0};           /*!< \brief Specific Heat Capacity at constant volume. */
-  su2double Mu{0.0};           /*!< \brief Laminar viscosity. */
-  su2double Mu_Turb{0.0};      /*!< \brief Eddy viscosity provided by a turbulence model. */
-  su2double dmudrho_T{0.0};    /*!< \brief Partial derivative of viscosity w.r.t. density. */
-  su2double dmudT_rho{0.0};    /*!< \brief Partial derivative of viscosity w.r.t. temperature. */
-  su2double Kt{0.0};           /*!< \brief Thermal conductivity. */
-  su2double dktdrho_T{0.0};    /*!< \brief Partial derivative of conductivity w.r.t. density. */
-  su2double dktdT_rho{0.0};    /*!< \brief Partial derivative of conductivity w.r.t. temperature. */
-  su2double mass_diffusivity{0.0};   /*!< \brief Mass Diffusivity */
+  su2double StaticEnergy{0.0};     /*!< \brief Internal Energy. */
+  su2double Entropy{0.0};          /*!< \brief Entropy. */
+  su2double Density{0.0};          /*!< \brief Density. */
+  su2double Pressure{0.0};         /*!< \brief Pressure. */
+  su2double SoundSpeed2{0.0};      /*!< \brief SpeedSound. */
+  su2double Temperature{0.0};      /*!< \brief Temperature. */
+  su2double dPdrho_e{0.0};         /*!< \brief DpDd_e. */
+  su2double dPde_rho{0.0};         /*!< \brief DpDe_d. */
+  su2double dTdrho_e{0.0};         /*!< \brief DTDd_e. */
+  su2double dTde_rho{0.0};         /*!< \brief DTDe_d. */
+  su2double dhdrho_P{0.0};         /*!< \brief DhDrho_p. */
+  su2double dhdP_rho{0.0};         /*!< \brief DhDp_rho. */
+  su2double dsdrho_P{0.0};         /*!< \brief DsDrho_p. */
+  su2double dsdP_rho{0.0};         /*!< \brief DsDp_rho. */
+  su2double Cp{0.0};               /*!< \brief Specific Heat Capacity at constant pressure. */
+  su2double Cv{0.0};               /*!< \brief Specific Heat Capacity at constant volume. */
+  su2double Mu{0.0};               /*!< \brief Laminar viscosity. */
+  su2double Mu_Turb{0.0};          /*!< \brief Eddy viscosity provided by a turbulence model. */
+  su2double dmudrho_T{0.0};        /*!< \brief Partial derivative of viscosity w.r.t. density. */
+  su2double dmudT_rho{0.0};        /*!< \brief Partial derivative of viscosity w.r.t. temperature. */
+  su2double Kt{0.0};               /*!< \brief Thermal conductivity. */
+  su2double dktdrho_T{0.0};        /*!< \brief Partial derivative of conductivity w.r.t. density. */
+  su2double dktdT_rho{0.0};        /*!< \brief Partial derivative of conductivity w.r.t. temperature. */
+  su2double mass_diffusivity{0.0}; /*!< \brief Mass Diffusivity */
 
   unique_ptr<CViscosityModel> LaminarViscosity;       /*!< \brief Laminar Viscosity Model */
   unique_ptr<CConductivityModel> ThermalConductivity; /*!< \brief Thermal Conductivity Model */
-  unique_ptr<CDiffusivityModel> MassDiffusivity;       /*!< \brief Mass Diffusivity Model */
+  unique_ptr<CDiffusivityModel> MassDiffusivity;      /*!< \brief Mass Diffusivity Model */
 
   /*!
    * \brief Instantiate the right type of viscosity model based on config.
@@ -147,7 +147,10 @@ class CFluidModel {
    * \brief Evaluate data manifold for flamelet or data-driven fluid problems.
    * \param[in] input - input data for manifold regression.
    */
-  virtual unsigned long EvaluateDataSet(const vector<su2double> &input_scalar, unsigned short lookup_type, vector<su2double> &output_refs) { return 0; }
+  virtual unsigned long EvaluateDataSet(const vector<su2double>& input_scalar, unsigned short lookup_type,
+                                        vector<su2double>& output_refs) {
+    return 0;
+  }
 
   /*!
    * \brief Get fluid dynamic viscosity.
@@ -331,7 +334,7 @@ class CFluidModel {
    * \brief Virtual member.
    * \param[in] T - Temperature value at the point.
    */
-  virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) { }
+  virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) {}
 
   /*!
    * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity.
@@ -344,6 +347,12 @@ class CFluidModel {
    */
   virtual unsigned long GetExtrapolation() const { return 0; }
 
+  /*!
+   * \brief Preferential diffusion as relevant phenomenon in flamelet simulations.
+   * \return Inclusion of preferential diffusion model.
+   */
+  virtual bool GetPreferentialDiffusion() const { return false; }
+
   /*!
    * \brief Get number of Newton solver iterations.
    * \return Newton solver iteration count at termination.

SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp

@@ -37,7 +37,6 @@
  */
 class CSpeciesFlameletSolver final : public CSpeciesSolver {
  private:
-  vector<su2activematrix> conjugate_var; /*!< \brief CHT variables for each boundary and vertex. */
   bool include_mixture_fraction = false; /*!< \brief include mixture fraction as a controlling variable. */
   /*!
    * \brief Compute the preconditioner for low-Mach flows.
@@ -63,8 +62,8 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \param[in] val_enth_out - pointer to output enthalpy variable.
    * \param[out] Converged - 0 if Newton solver converged, 1 if not.
    */
-  unsigned long GetEnthFromTemp(CFluidModel* fluid_model, su2double const val_temp,
-                                const su2double* scalar_solution, su2double* val_enth_out);
+  unsigned long GetEnthFromTemp(CFluidModel* fluid_model, su2double const val_temp, const su2double* scalar_solution,
+                                su2double* val_enth_out);
 
   /*!
    * \brief Find maximum progress variable value within the manifold for the current solution.
@@ -97,13 +96,23 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
   unsigned long SetScalarLookUps(const CConfig* config, CFluidModel* fluid_model_local, unsigned long iPoint,
                                  const vector<su2double>& scalars);
 
+  /*!
+   * \brief Retrieve the preferential diffusion scalar values from manifold.
+   * \param[in] config - definition of particular problem.
+   * \param[in] fluid_model_local - pointer to flamelet fluid model.
+   * \param[in] iPoint - node ID.
+   * \param[in] scalars - local scalar solution.
+   * \return - within manifold bounds (0) or outside manifold bounds (1).
+   */
+  unsigned long SetPreferentialDiffusionScalars(const CConfig* config, CFluidModel* fluid_model_local,
+                                                unsigned long iPoint, const vector<su2double>& scalars);
+
  public:
   /*!
-   * \brief Constructor.
+   * \brief Define a Flamelet Generated Manifold species solver.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Index of the mesh in multigrid computations.
-   * \param[in] FluidModel
    */
   CSpeciesFlameletSolver(CGeometry* geometry, CConfig* config, unsigned short iMesh);
 
@@ -175,4 +184,16 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    */
   void BC_ConjugateHeat_Interface(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, CConfig* config,
                                   unsigned short val_marker) override;
+
+  /*!
+   * \brief Compute the fluxes due to viscous and preferential diffusion effects of the flamelet species at a particular edge.
+   * \param[in] iEdge - Edge for which we want to compute the flux
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \note Calls a generic implementation after defining a SolverSpecificNumerics object.
+   */
+  void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, CNumerics* numerics,
+                        CConfig* config) final;
 };