[WIP] Assumed shape beta pdf

Common/include/option_structure.hpp

@@ -1525,14 +1525,23 @@ static const MapType<std::string, SPECIES_MODEL> Species_Model_Map = {
 };
 
 /*!
- * \brief Progress variable and enthalpy are the first and second entries in the lookup table.
- * \note The order matters.
+ * \brief Names of the flamelet controlling variables in CONTROLLING_VARIABLE_NAMES; the name sets the role.
+ */
+namespace FLAMELET_CV_NAME {
+  constexpr const char* PROGVAR = "ProgressVariable";
+  constexpr const char* ENTHALPY = "EnthalpyTot";
+  constexpr const char* MIXFRAC = "MixtureFraction";
+  constexpr const char* PROGVAR_VAR = "ProgressVariableVariance";
+  constexpr const char* MIXFRAC_VAR = "MixtureFractionVariance";
+}
+
+/*!
+ * \brief Manifold variables giving the range of the progress variable, the support of its PDF.
  */
-enum FLAMELET_SCALAR_VARIABLES {
-  I_PROGVAR,
-  I_ENTH,
-  I_MIXFRAC,
-};
+namespace FLAMELET_PROGVAR_BOUNDS {
+  constexpr const char* MIN = "ProgressVariableMin";
+  constexpr const char* MAX = "ProgressVariableMax";
+}
 
 /*!
  * \brief the source terms for the flamelet method. At the moment only progress variable
@@ -1549,6 +1558,7 @@ enum FLAMELET_LOOKUP_OPS {
   PREFDIF,  /*!< \brief Preferential diffusion scalars. */
   SOURCES,  /*!< \brief Scalar source terms (controlling variables, passive species).*/
   LOOKUP,   /*!< \brief Passive look-up variables specified in config. */
+  PROGVAR_BOUNDS, /*!< \brief Minimum and maximum of the progress variable. */
 };
 
 /*!
@@ -1607,6 +1617,12 @@ struct FluidFlamelet_ParsedOptions {
   unsigned short n_user_scalars = 0;  /*!< \brief Number of user defined (auxiliary) scalar transport equations. */
   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. */
+  int idx_progvar = -1;               /*!< \brief Scalar index of the progress variable, -1 if absent. */
+  int idx_enth = -1;                  /*!< \brief Scalar index of the total enthalpy. */
+  int idx_mixfrac = -1;               /*!< \brief Scalar index of the mixture fraction, -1 if absent. */
+  int idx_progvar_var = -1;           /*!< \brief Scalar index of the progress variable variance, -1 if absent. */
+  int idx_mixfrac_var = -1;           /*!< \brief Scalar index of the mixture fraction variance, -1 if absent. */
+  su2double variance_dissipation_coeff = 2.0; /*!< \brief C_d in the variance dissipation C_d rho (eps/k) var. */
 
   std::string *controlling_variable_names; /*!< \brief Names of the independent, transported scalars. */
   std::string* cv_source_names;            /*!< \brief Names of the source terms of the independent, transported scalars. */

Common/src/CConfig.cpp

@@ -1528,6 +1528,9 @@ void CConfig::SetConfig_Options() {
   /*!\brief FLAME_INIT_IGNITION \n DESCRIPTION: Ignition temperature for the flame initialization \ingroup Config*/
   addDoubleOption("FLAME_INIT_IGNITION", flamelet_ParsedOptions.Flame_T_ignition, 5000.0);
 
+  /*!\brief FLAMELET_VARIANCE_DISSIPATION_COEFF \n DESCRIPTION: Coefficient C_d of the scalar variance dissipation C_d rho (eps/k) var \ingroup Config*/
+  addDoubleOption("FLAMELET_VARIANCE_DISSIPATION_COEFF", flamelet_ParsedOptions.variance_dissipation_coeff, 2.0);
+
   /*--- Options related to mass diffusivity and thereby the species solver. ---*/
 
   /*!\brief DIFFUSIVITY_MODEL\n DESCRIPTION: mass diffusivity model \n DEFAULT constant disffusivity \ingroup Config*/
@@ -6109,6 +6112,51 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     /*--- We can have additional user defined transported scalars ---*/
     flamelet_ParsedOptions.n_scalars = flamelet_ParsedOptions.n_control_vars + flamelet_ParsedOptions.n_user_scalars;
 
+    /*--- The name of each controlling variable sets its role. ---*/
+    auto& opts = flamelet_ParsedOptions;
+    for (int iCV = 0; iCV < opts.n_control_vars; iCV++) {
+      const string& name = opts.controlling_variable_names[iCV];
+      int* idx = nullptr;
+      if (name == FLAMELET_CV_NAME::PROGVAR) idx = &opts.idx_progvar;
+      else if (name == FLAMELET_CV_NAME::ENTHALPY) idx = &opts.idx_enth;
+      else if (name == FLAMELET_CV_NAME::MIXFRAC) idx = &opts.idx_mixfrac;
+      else if (name == FLAMELET_CV_NAME::PROGVAR_VAR) idx = &opts.idx_progvar_var;
+      else if (name == FLAMELET_CV_NAME::MIXFRAC_VAR) idx = &opts.idx_mixfrac_var;
+      else
+        SU2_MPI::Error("Unknown controlling variable " + name + ", expected " + FLAMELET_CV_NAME::PROGVAR + ", " +
+                       FLAMELET_CV_NAME::ENTHALPY + ", " + FLAMELET_CV_NAME::MIXFRAC + ", " +
+                       FLAMELET_CV_NAME::PROGVAR_VAR + " or " + FLAMELET_CV_NAME::MIXFRAC_VAR + ".",
+                       CURRENT_FUNCTION);
+      if (*idx >= 0)
+        SU2_MPI::Error("Controlling variable " + name + " is listed twice.",
+                       CURRENT_FUNCTION);
+      *idx = iCV;
+    }
+    if (opts.idx_enth < 0)
+      SU2_MPI::Error(string("CONTROLLING_VARIABLE_NAMES must include ") + FLAMELET_CV_NAME::ENTHALPY + ".",
+                     CURRENT_FUNCTION);
+    if (opts.idx_progvar < 0 && opts.idx_mixfrac < 0)
+      SU2_MPI::Error(string("CONTROLLING_VARIABLE_NAMES must include ") + FLAMELET_CV_NAME::PROGVAR + " or " +
+                     FLAMELET_CV_NAME::MIXFRAC + ".", CURRENT_FUNCTION);
+
+    /*--- The thickened flame correction measures the progress variable gradient. ---*/
+    if (opts.idx_progvar < 0) opts.thickenedflame_correction = false;
+
+    /*--- A variance needs its mean and the turbulence time scale omega of the SST model. ---*/
+    if ((opts.idx_progvar_var >= 0) && (opts.idx_progvar < 0))
+      SU2_MPI::Error(string(FLAMELET_CV_NAME::PROGVAR_VAR) + " requires " + FLAMELET_CV_NAME::PROGVAR + ".",
+                     CURRENT_FUNCTION);
+    if ((opts.idx_mixfrac_var >= 0) && (opts.idx_mixfrac < 0))
+      SU2_MPI::Error(string(FLAMELET_CV_NAME::MIXFRAC_VAR) + " requires " + FLAMELET_CV_NAME::MIXFRAC + ".",
+                     CURRENT_FUNCTION);
+    if ((opts.idx_progvar_var >= 0) || (opts.idx_mixfrac_var >= 0)) {
+      if (Kind_Turb_Model != TURB_MODEL::SST)
+        SU2_MPI::Error("Scalar variance transport requires KIND_TURB_MODEL= SST.", CURRENT_FUNCTION);
+      if (opts.preferential_diffusion)
+        SU2_MPI::Error("PREFERENTIAL_DIFFUSION is not available with scalar variance transport.", CURRENT_FUNCTION);
+      opts.thickenedflame_correction = false;
+    }
+
     /*--- Check that spark ignition has required parameters defined ---*/
     if (flamelet_ParsedOptions.ignition_method == FLAMELET_INIT_TYPE::SPARK) {
       /*--- Check if SPARK_INIT was explicitly set in config file ---*/

SU2_CFD/include/fluid/CFluidFlamelet.hpp

@@ -48,8 +48,7 @@ class CFluidFlamelet final : public CFluidModel {
   su2double beta_progvar, beta_enth_thermal, beta_enth, beta_mixfrac;
   int rank;
 
-  bool include_mixture_fraction = false, /*!< \brief include mixture fraction in controlling variables. */
-       preferential_diffusion = false;     /*!< \brief use preferential diffusion physics. */
+  bool preferential_diffusion = 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. */
@@ -68,7 +67,12 @@ class CFluidFlamelet final : public CFluidModel {
   vector<unsigned long> LUT_idx_TD,
                         LUT_idx_Sources,
                         LUT_idx_LookUp,
-                        LUT_idx_PD;
+                        LUT_idx_PD,
+                        LUT_idx_ProgVarBounds;
+
+  vector<su2double> manifold_inputs; /*!< \brief Controlling variables with the variances normalized. */
+  su2double limit_progvar_var = 0,   /*!< \brief Upper limit of the progress variable variance. */
+      limit_mixfrac_var = 0;         /*!< \brief Upper limit of the mixture fraction variance. */
 
   /*--- Class variables for the multi-layer perceptron method ---*/
 #ifdef USE_MLPCPP
@@ -78,23 +82,41 @@ class CFluidFlamelet final : public CFluidModel {
   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;
+  MLPToolbox::CIOMap* iomap_ProgVarBounds; /*!< \brief Input-output map for the progress variable range. */
 #endif
 
   vector<su2double> scalars_vector;
 
   vector<string> varnames_TD, /*!< \brief Lookup names for thermodynamic state variables. */
       varnames_Sources,       /*!< \brief Lookup names for source terms. */
       varnames_LookUp,        /*!< \brief Lookup names for passive look-up terms. */
-      varnames_PD;            /*!< \brief Lookup names for preferential diffusion scalars. */
+      varnames_PD,            /*!< \brief Lookup names for preferential diffusion scalars. */
+      varnames_ProgVarBounds; /*!< \brief Lookup names for the minimum and maximum progress variable. */
 
   vector<su2double> val_vars_TD, /*!< \brief References to thermodynamic state variables. */
       val_vars_Sources,          /*!< \brief References to source terms. */
       val_vars_LookUp,           /*!< \brief References passive look-up terms. */
-      val_vars_PD;               /*!< \brief References to preferential diffusion scalars. */
+      val_vars_PD,               /*!< \brief References to preferential diffusion scalars. */
+      val_vars_ProgVarBounds;    /*!< \brief References to the minimum and maximum progress variable. */
 
   void PreprocessLookUp(CConfig* config);
 
+  /*!
+   * \brief Interpolate the manifold, without normalizing the inputs.
+   * \param[in] lookup_type - Look-up operation to be performed (FLAMELET_LOOKUP_OPS).
+   * \param[in] inputs - Manifold inputs in the order of the controlling variables.
+   * \param[out] output_refs - Interpolated output variables.
+   * \returns Whether the query is inside the manifold.
+   */
+  bool LookUpManifold(unsigned short lookup_type, const vector<su2double>& inputs, vector<su2double>& output_refs);
+
+  /*!
+   * \brief Set manifold_inputs to the controlling variables, each variance divided by its upper limit.
+   * \param[in] scalars - Transported scalars.
+   * \returns Whether the look-up of the progress variable range is inside the manifold.
+   */
+  bool SetManifoldInputs(const su2double* scalars);
+
   /*! \brief
    * Returns true if the string is null or zero (ignores case).
    */
@@ -135,6 +157,14 @@ class CFluidFlamelet final : public CFluidModel {
   unsigned long EvaluateDataSet(const vector<su2double>& input_scalar, unsigned short lookup_type,
                                        vector<su2double>& output_refs);
 
+  /*!
+   * \brief Largest value a scalar variance can take, (mean - min)(max - mean) for its mean.
+   * \param[in] scalars - Transported scalars.
+   * \param[in] i_var - Index of the variance.
+   * \returns Upper limit of the variance.
+   */
+  su2double GetScalarVarianceLimit(const su2double* scalars, unsigned short i_var) override;
+
   /*!
    * \brief Check for out-of-bounds condition for data set interpolation.
    * \return - within bounds (0) or out of bounds (1).

SU2_CFD/include/fluid/CFluidModel.hpp

@@ -161,6 +161,14 @@ class CFluidModel {
     return 0;
   }
 
+  /*!
+   * \brief Largest value a scalar variance can take, (mean - min)(max - mean) for its mean.
+   * \param[in] scalars - Transported scalars.
+   * \param[in] i_var - Index of the variance.
+   * \returns Upper limit of the variance.
+   */
+  virtual su2double GetScalarVarianceLimit(const su2double* scalars, unsigned short i_var) { return 1.0e15; }
+
   /*!
    * \brief Get fluid dynamic viscosity.
    */

SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp

@@ -50,7 +50,9 @@ class CScalarFlux_Flamelet final
   explicit CScalarFlux_Flamelet(const CConfig& config)
       : Base(config),
         preferentialDiffusion(config.GetFlameletParsedOptions().preferential_diffusion),
-        nControlVars(config.GetFlameletParsedOptions().n_control_vars) {}
+        nControlVars(config.GetFlameletParsedOptions().n_control_vars),
+        idxProgVar(config.GetFlameletParsedOptions().idx_progvar),
+        idxEnth(config.GetFlameletParsedOptions().idx_enth) {}
 
   /*!
    * \brief Preferential diffusion, two terms with the shape of the ordinary diffusion but of
@@ -114,42 +116,36 @@ class CScalarFlux_Flamelet final
     }
 
     /*--- Thermal term, on the enthalpy equation alone, driven by the temperature gradient. ---*/
-    if (nControlVars <= I_ENTH) return;
-
     const Double T_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.Temperature());
     const Double T_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.Temperature());
 
     const auto gradT_i = gatherVariables<nDim>(iPoint, side_i.flowNodes->GetGradient_Primitive(), idx.Temperature());
     const auto gradT_j = gatherVariables<nDim>(jPoint, side_j.flowNodes->GetGradient_Primitive(), idx.Temperature());
 
     const Double Dth_i = gatherVariables(iPoint, side_i.scalarNodes.GetAuxVar(), I_BETA_ENTH_THERMAL) *
-                         gatherVariables(iPoint, side_i.scalarNodes.GetDiffusivity(), I_ENTH);
+                         gatherVariables(iPoint, side_i.scalarNodes.GetDiffusivity(), idxEnth);
     const Double Dth_j = gatherVariables(jPoint, side_j.scalarNodes.GetAuxVar(), I_BETA_ENTH_THERMAL) *
-                         gatherVariables(jPoint, side_j.scalarNodes.GetDiffusivity(), I_ENTH);
+                         gatherVariables(jPoint, side_j.scalarNodes.GetDiffusivity(), idxEnth);
     const Double Dth = 0.5 * (rho.i * Dth_i + rho.j * Dth_j) + diffTurb;
 
     const Double projGradT = projectedGradient(opt, gradT_i, gradT_j, T_i, T_j, normal, vector_ij, dist2_ij);
 
-    res.flux_i(I_ENTH) -= Dth * projGradT;
-    if (!opt.oneSided) res.flux_j(I_ENTH) += Dth * projGradT;
+    res.flux_i(idxEnth) -= Dth * projGradT;
+    if (!opt.oneSided) res.flux_j(idxEnth) += Dth * projGradT;
   }
 
  private:
   const bool preferentialDiffusion;
   const unsigned short nControlVars;
+  const int idxProgVar, idxEnth;
 
   /*!
    * \brief Auxiliary variable holding the beta scalar of a controlling variable.
    */
-  static FORCEINLINE unsigned short betaIndex(unsigned short iScalar) {
-    switch (iScalar) {
-      case I_PROGVAR:
-        return I_BETA_PROGVAR;
-      case I_ENTH:
-        return I_BETA_ENTH;
-      default:
-        return I_BETA_MIXFRAC;
-    }
+  FORCEINLINE unsigned short betaIndex(int iScalar) const {
+    if (iScalar == idxProgVar) return I_BETA_PROGVAR;
+    if (iScalar == idxEnth) return I_BETA_ENTH;
+    return I_BETA_MIXFRAC;
   }
 
   /*!

SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp

@@ -41,7 +41,10 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
   const su2double default_flame_thickness{1.0};
   su2double global_flame_thickness;
   bool calc_flame_thickness{false};
-  bool include_mixture_fraction = false; /*!< \brief include mixture fraction as a controlling variable. */
+  int idx_progvar = -1; /*!< \brief Scalar index of the progress variable, -1 if absent. */
+  int idx_enth = -1;    /*!< \brief Scalar index of the total enthalpy. */
+  int idx_flame_front = -1; /*!< \brief Scalar set by FLAME_FRONT: progress variable, else mixture fraction. */
+  vector<pair<int, int>> variance_sources; /*!< \brief Variance and mean indices of variances with a chemical source. */
   /*!< \brief Number of points outside the manifold domain, shared across OpenMP threads so it can be
    *          accumulated atomically and reduced by the master thread alone in Preprocessing. */
   unsigned long n_not_in_domain_local = 0;
@@ -51,8 +54,27 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
+  /*!
+   * \brief Limit the scalar variances to the maximum of a two-point distribution on the range of their mean.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   */
+  void Postprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config,
+                      unsigned short iMesh) override;
+
   void SetPreconditioner(CGeometry* geometry, CSolver** solver_container, CConfig* config);
 
+  /*!
+   * \brief Add the turbulent production 2 mu_t/Sc_t |grad(mean)|^2 and dissipation C_d rho beta* omega var
+   *        of a scalar variance to the residual, the dissipation implicitly.
+   * \param[in] i_var - Scalar index of the variance.
+   * \param[in] i_mean - Scalar index of its mean.
+   */
+  void ScalarVarianceResidual(CGeometry* geometry, CSolver** solver_container, const CConfig* config, int i_var,
+                              int i_mean);
+
   /*!
    * \brief Generic implementation of the isothermal wall also covering CHT cases,
    * for which the wall temperature is given by GetConjugateHeatVariable.