Improved combustion solver robustness for coarse grids

Common/include/option_structure.hpp

@@ -1448,6 +1448,8 @@ struct FluidFlamelet_ParsedOptions {
   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.*/
+
+  su2double flame_lengthscale{1e-3}; /*!< \brief Lengthscale above which chemical source terms are dampened.*/
 };
 
 /*!

Common/src/CConfig.cpp

@@ -1433,6 +1433,9 @@ void CConfig::SetConfig_Options() {
   /*!\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", flamelet_ParsedOptions.nspark, flamelet_ParsedOptions.spark_reaction_rates);
 
+  /*!\brief FLAME_LENGTHSCALE \n DESCRIPTION: Lengthscale above which species source terms are dampened. \ingroup Config*/
+  addDoubleOption("FLAME_LENGTHSCALE", flamelet_ParsedOptions.flame_lengthscale, 1e-3);
+
   /*--- Options related to mass diffusivity and thereby the species solver. ---*/
 
   /*!\brief DIFFUSIVITY_MODEL\n DESCRIPTION: mass diffusivity model \n DEFAULT constant disffusivity \ingroup Config*/
@@ -5797,6 +5800,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Number of initial species incompatible with number of controlling variables and user scalars.", CURRENT_FUNCTION);
     /*--- We can have additional user defined transported scalars ---*/
     flamelet_ParsedOptions.n_scalars = flamelet_ParsedOptions.n_control_vars + flamelet_ParsedOptions.n_user_scalars;
+
+    if (flamelet_ParsedOptions.flame_lengthscale <= 0.0)
+      SU2_MPI::Error("Flame length scale value should be positive.", CURRENT_FUNCTION);
+
   }
 
   if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) {

SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp

@@ -213,7 +213,6 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver**
 
     for (unsigned long i_mesh = 0; i_mesh <= config->GetnMGLevels(); i_mesh++) {
       fluid_model_local = solver_container[i_mesh][FLOW_SOL]->GetFluidModel();
-      prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init);
       for (auto iVar = 0u; iVar < nVar; iVar++) scalar_init[iVar] = config->GetSpecies_Init()[iVar];
 
       /*--- Set enthalpy based on initial temperature and scalars. ---*/
@@ -227,6 +226,7 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver**
 
         if (flame_front_ignition) {
 
+          prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init);
           /*--- Determine if point is above or below the plane, assuming the normal
             is pointing towards the burned region. ---*/
           point_loc = 0.0;
@@ -382,12 +382,16 @@ void CSpeciesFlameletSolver::SetPreconditioner(CGeometry* geometry, CSolver** so
 
 void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solver_container,
                                              CNumerics** numerics_container, CConfig* config, unsigned short iMesh) {
+
+   const su2double flame_lengthscale = config->GetFlameletParsedOptions().flame_lengthscale;  
+   const su2double flame_vol = pow(flame_lengthscale, nDim); 
 
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto i_point = 0u; i_point < nPointDomain; i_point++) {
+    const su2double fac = min(flame_vol/ geometry->nodes->GetVolume(i_point), 1.0);
     /*--- Add source terms from the lookup table directly to the residual. ---*/
     for (auto i_var = 0; i_var < nVar; i_var++) {
-      LinSysRes(i_point, i_var) -= nodes->GetScalarSources(i_point)[i_var] * geometry->nodes->GetVolume(i_point);
+      LinSysRes(i_point, i_var) -= fac*nodes->GetScalarSources(i_point)[i_var] * geometry->nodes->GetVolume(i_point);
     }
   }
   END_SU2_OMP_FOR