From 11a88fb506dcd8561817a7ea79cd461084eda637 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Tue, 17 Feb 2026 14:12:49 +0100
Subject: [PATCH 01/23] Merge branch 'develop' of
 https://github.com/su2code/SU2 into feature_flamelet

---
 Common/include/option_structure.hpp            | 2 ++
 Common/src/CConfig.cpp                         | 7 +++++++
 SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp | 8 ++++++--
 3 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 690350fae05..e9453b6dd0e 100644
--- a/Common/include/option_structure.hpp
+++ b/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.*/
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 95f5754ef34..5b9582bc776 100644
--- a/Common/src/CConfig.cpp
+++ b/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()) {
diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 622a1a6b068..e72eab09a3e 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/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

From 6ba1ffb2538b3edcf2446e69491e50540e451432 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Tue, 17 Feb 2026 14:27:04 +0100
Subject: [PATCH 02/23] fix trailing whitespace

---
 Common/src/CConfig.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 5b9582bc776..ff6bf6aee42 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -5801,9 +5801,9 @@ 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;
 
-    if (flamelet_ParsedOptions.flame_lengthscale <= 0.0) 
+    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()) {

From 2df59a8e488ac4d2b8643081b649c4db7cff1535 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Fri, 20 Feb 2026 10:53:40 +0100
Subject: [PATCH 03/23] Added function for thickened flame correction factor

---
 SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
index df61813a3e2..26f44286691 100644
--- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
+++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
@@ -81,10 +81,11 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \param[in] iPoint - node ID.
    * \param[in] scalars - local scalar solution.
    * \param[in] table_source_names - variable names of scalar source terms.
+   * \param[in] F - flame thickness correction factor.
    * \return - within manifold bounds (0) or outside manifold bounds (1).
    */
   unsigned long SetScalarSources(const CConfig* config, CFluidModel* fluid_model_local, unsigned long iPoint,
-                                 const vector<su2double>& scalars);
+                                 const vector<su2double>& scalars, const su2double F=1.0);
 
   /*!
    * \brief Retrieve passive look-up data from manifold.
@@ -105,6 +106,14 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    */
   unsigned long SetPreferentialDiffusionScalars(CFluidModel* fluid_model_local,
                                                 unsigned long iPoint, const vector<su2double>& scalars);
+  
+  /*!
+   * \brief Calculate correction factor for flame propagation on coarse grids.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] iPoint - node ID.
+   * \return - flame thickness correction factor.
+   */                                              
+  su2double ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint);   
 
  public:
   /*!

From 3af071244848301325d3bd7a59043feafbe1d759 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Fri, 20 Feb 2026 10:55:02 +0100
Subject: [PATCH 04/23] Diffusivity is multiplied and source terms are divided
 by correction factor according to the model by Butler and O'Rouke

---
 .../src/solvers/CSpeciesFlameletSolver.cpp    | 30 ++++++++++++-------
 1 file changed, 19 insertions(+), 11 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index e72eab09a3e..464336d0a56 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -91,10 +91,14 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
   for (auto i_point = 0u; i_point < nPoint; i_point++) {
     CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel();
     su2double* scalars = nodes->GetSolution(i_point);
+
+    /*--- Calculate correction factor for flame propagation on coarse grids. ---*/
+    const su2double F = ThickenedFlameCorrection(geometry, i_point);
+
     for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar];
 
     /*--- Compute total source terms from the production and consumption. ---*/
-    unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector);
+    unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector, F);
 
     if (ignition) {
       /*--- Apply source terms within spark radius. ---*/
@@ -103,7 +107,7 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
       dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), flamelet_config_options.spark_init.data());
       if (dist_from_center < pow(spark_radius,2)) {
         for (auto iVar = 0u; iVar < nVar; iVar++)
-          nodes->SetScalarSource(i_point, iVar, nodes->GetScalarSources(i_point)[iVar] + flamelet_config_options.spark_reaction_rates[iVar]);
+          nodes->SetScalarSource(i_point, iVar,  nodes->GetScalarSources(i_point)[iVar] + flamelet_config_options.spark_reaction_rates[iVar] / F);
       }
     }
 
@@ -115,9 +119,9 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
     /*--- Set mass diffusivity based on thermodynamic state. ---*/
     auto T = flowNodes->GetTemperature(i_point);
     fluid_model_local->SetTDState_T(T, scalars);
-    /*--- set the diffusivity in the fluid model to the diffusivity obtained from the lookup table ---*/
+    /*--- set the diffusivity in the fluid model to the diffusivity obtained from the lookup table, multiplied by flame thickness correction factor ---*/
     for (auto i_scalar = 0u; i_scalar < nVar; ++i_scalar) {
-      nodes->SetDiffusivity(i_point, fluid_model_local->GetMassDiffusivity(i_scalar), i_scalar);
+      nodes->SetDiffusivity(i_point, F * (fluid_model_local->GetMassDiffusivity(i_scalar)), i_scalar);
     }
 
     /*--- Obtain preferential diffusion scalar values. ---*/
@@ -383,15 +387,11 @@ 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) -= fac*nodes->GetScalarSources(i_point)[i_var] * geometry->nodes->GetVolume(i_point);
+      LinSysRes(i_point, i_var) -= nodes->GetScalarSources(i_point)[i_var] * geometry->nodes->GetVolume(i_point);
     }
   }
   END_SU2_OMP_FOR
@@ -519,7 +519,7 @@ void CSpeciesFlameletSolver::BC_ConjugateHeat_Interface(CGeometry* geometry, CSo
 }
 
 unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CFluidModel* fluid_model_local,
-                                                       unsigned long iPoint, const vector<su2double>& scalars) {
+                                                       unsigned long iPoint, const vector<su2double>& scalars, const su2double F) {
   /*--- Compute total source terms from the production and consumption. ---*/
 
   vector<su2double> table_sources(flamelet_config_options.n_control_vars + 2 * flamelet_config_options.n_user_scalars);
@@ -541,8 +541,9 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF
     su2double source_cons = table_sources[flamelet_config_options.n_control_vars + 2 * i_aux + 1];
     source_scalar[flamelet_config_options.n_control_vars + i_aux] = source_prod + source_cons * y_aux;
   }
+  /*--- Source term is divided by flame thickness correction factor to improve stability on coarse grids. ---*/
   for (auto i_scalar = 0u; i_scalar < nVar; i_scalar++)
-    nodes->SetScalarSource(iPoint, i_scalar, source_scalar[i_scalar]);
+    nodes->SetScalarSource(iPoint, i_scalar,  source_scalar[i_scalar] / F);
   return misses;
 }
 
@@ -807,3 +808,10 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo
   su2double pv_burnt = scalars[I_PROGVAR] - delta;
   return pv_burnt;
 }
+
+
+su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) {
+  su2double flame_vol = pow(flamelet_config_options.flame_lengthscale,nDim);
+  su2double F = max(1.0, geometry->nodes->GetVolume(iPoint) / flame_vol);
+  return F;
+}
\ No newline at end of file

From e361ca705b0c3ddba3d8d2589ecb6025855fbb55 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Fri, 20 Feb 2026 13:17:38 +0100
Subject: [PATCH 05/23] Apply source correction only for steady problems

---
 SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 464336d0a56..865fe33f8d4 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -81,7 +81,15 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
     auto spark_init = flamelet_config_options.spark_init;
     spark_iter_start = ceil(spark_init[4]);
     spark_duration = ceil(spark_init[5]);
-    unsigned long iter = config->GetMultizone_Problem() ? config->GetOuterIter() : config->GetInnerIter();
+    unsigned long iter;
+    if (config->GetMultizone_Problem()) {
+      iter = config->GetOuterIter();
+    } else if (config->GetTime_Domain()) {
+      iter = config->GetTimeIter();
+    } else {
+      iter = config->GetInnerIter();
+    }
+    //unsigned long iter = (config->GetMultizone_Problem() || config->GetTime_Domain()) ? config->GetOuterIter() : config->GetInnerIter();
     ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)));
   }
 
@@ -98,7 +106,10 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
     for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar];
 
     /*--- Compute total source terms from the production and consumption. ---*/
-    unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector, F);
+
+    /*--- Only apply thickened flame correction factor to sources for steady problems. ---*/
+    su2double F_source = config->GetTime_Domain() ? 1.0 : 1.0 / F;
+    unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector, F_source);
 
     if (ignition) {
       /*--- Apply source terms within spark radius. ---*/

From 0de7e0dd8b379f5d1be2d7b5f130bf72901332cb Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Fri, 20 Feb 2026 13:43:49 +0100
Subject: [PATCH 06/23] bug fix

---
 SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 865fe33f8d4..e2548c19158 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -72,6 +72,7 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
                                            unsigned short RunTime_EqSystem, bool Output) {
   unsigned long n_not_in_domain_local = 0, n_not_in_domain_global = 0;
   vector<su2double> scalars_vector(nVar);
+  const su2double F_default{1.0};
   unsigned long spark_iter_start, spark_duration;
   bool ignition = false;
   auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
@@ -89,7 +90,6 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
     } else {
       iter = config->GetInnerIter();
     }
-    //unsigned long iter = (config->GetMultizone_Problem() || config->GetTime_Domain()) ? config->GetOuterIter() : config->GetInnerIter();
     ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)));
   }
 
@@ -101,14 +101,14 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
     su2double* scalars = nodes->GetSolution(i_point);
 
     /*--- Calculate correction factor for flame propagation on coarse grids. ---*/
-    const su2double F = ThickenedFlameCorrection(geometry, i_point);
+    su2double F = ThickenedFlameCorrection(geometry, i_point);
 
     for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar];
 
     /*--- Compute total source terms from the production and consumption. ---*/
 
     /*--- Only apply thickened flame correction factor to sources for steady problems. ---*/
-    su2double F_source = config->GetTime_Domain() ? 1.0 : 1.0 / F;
+    su2double F_source = config->GetTime_Domain() ? F_default : 1.0 / F;
     unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector, F_source);
 
     if (ignition) {
@@ -554,7 +554,7 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF
   }
   /*--- Source term is divided by flame thickness correction factor to improve stability on coarse grids. ---*/
   for (auto i_scalar = 0u; i_scalar < nVar; i_scalar++)
-    nodes->SetScalarSource(iPoint, i_scalar,  source_scalar[i_scalar] / F);
+    nodes->SetScalarSource(iPoint, i_scalar,  F*source_scalar[i_scalar]);
   return misses;
 }
 

From ff6faea69dc11a19deec4f0f534b15536f84ab58 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Thu, 19 Mar 2026 11:04:37 +0100
Subject: [PATCH 07/23] NULL MLP no longer needed

---
 .../laminar_premixed_h2_flame_cfd.cfg                           | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg b/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg
index 93328224dc0..974505d809b 100644
--- a/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg
+++ b/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg
@@ -31,7 +31,7 @@ CONDUCTIVITY_MODEL= FLAMELET
 DIFFUSIVITY_MODEL= FLAMELET
 KIND_SCALAR_MODEL= FLAMELET
 INTERPOLATION_METHOD= MLP
-FILENAMES_INTERPOLATOR= (MLP_TD.mlp, MLP_PD.mlp, MLP_PPV.mlp, MLP_null.mlp)
+FILENAMES_INTERPOLATOR= (MLP_TD.mlp, MLP_PD.mlp, MLP_PPV.mlp)
 PREFERENTIAL_DIFFUSION= YES
 
 % -------------------- SCALAR TRANSPORT ---------------------------------------%

From 5a9d2a0512a79b6d633e5d7f79f227873159e1d9 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Thu, 19 Mar 2026 11:08:13 +0100
Subject: [PATCH 08/23] Added general flame thickness calculation method

---
 .../src/solvers/CSpeciesFlameletSolver.cpp    | 100 ++++++++++++++----
 1 file changed, 82 insertions(+), 18 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index e2548c19158..8161c051486 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -40,6 +40,14 @@ CSpeciesFlameletSolver::CSpeciesFlameletSolver(CGeometry* geometry, CConfig* con
   /*--- Retrieve options from config. ---*/
   flamelet_config_options = config->GetFlameletParsedOptions();
 
+  if (flamelet_config_options.flame_lengthscale != flamelet_config_options.default_lengthscale) {
+    global_flame_thickness = flamelet_config_options.flame_lengthscale;
+    calc_flame_thickness = false;
+  } else {
+    global_flame_thickness = default_flame_thickness;
+    calc_flame_thickness = true;
+  }
+
   /*--- Dimension of the problem. ---*/
   nVar = flamelet_config_options.n_scalars;
   include_mixture_fraction = (flamelet_config_options.n_control_vars == 3);
@@ -65,6 +73,14 @@ CSpeciesFlameletSolver::CSpeciesFlameletSolver(CGeometry* geometry, CConfig* con
 
   /*--- Add the solver name. ---*/
   SolverName = "FLAMELET";
+
+  if (rank==MASTER_NODE) {
+    if (calc_flame_thickness) {
+      cout << "Lengthscale for thickened flame model is automatically determined." << endl;
+    } else {
+      cout << "Thickened flame model is applied in cells larger than " << global_flame_thickness*1000 << " mm" << endl;
+    }
+  }
 }
 
 void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config,
@@ -72,43 +88,46 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
                                            unsigned short RunTime_EqSystem, bool Output) {
   unsigned long n_not_in_domain_local = 0, n_not_in_domain_global = 0;
   vector<su2double> scalars_vector(nVar);
-  const su2double F_default{1.0};
+  
   unsigned long spark_iter_start, spark_duration;
   bool ignition = false;
   auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
   /*--- Retrieve spark ignition parameters for spark-type ignition. ---*/
+  unsigned long iter;
+  if (config->GetMultizone_Problem()) {
+    iter = config->GetOuterIter();
+  } else if (config->GetTime_Domain()) {
+    iter = config->GetTimeIter();
+  } else {
+    iter = config->GetInnerIter();
+  }
   if ((flamelet_config_options.ignition_method == FLAMELET_INIT_TYPE::SPARK)) {
     auto spark_init = flamelet_config_options.spark_init;
     spark_iter_start = ceil(spark_init[4]);
     spark_duration = ceil(spark_init[5]);
-    unsigned long iter;
-    if (config->GetMultizone_Problem()) {
-      iter = config->GetOuterIter();
-    } else if (config->GetTime_Domain()) {
-      iter = config->GetTimeIter();
-    } else {
-      iter = config->GetInnerIter();
-    }
+    
     ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)));
   }
-
   SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);)
 
+  su2double test_thickness = GetOverallFlameThickness(geometry, solver_container);
+  if (test_thickness < global_flame_thickness) {
+    global_flame_thickness = test_thickness;
+  } 
+  if (rank==MASTER_NODE) cout << global_flame_thickness << endl;
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto i_point = 0u; i_point < nPoint; i_point++) {
     CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel();
     su2double* scalars = nodes->GetSolution(i_point);
 
     /*--- Calculate correction factor for flame propagation on coarse grids. ---*/
-    su2double F = ThickenedFlameCorrection(geometry, i_point);
+    const su2double F = ThickenedFlameCorrection(geometry, i_point);
+    const su2double F_source = 1.0 / F;
 
     for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar];
 
-    /*--- Compute total source terms from the production and consumption. ---*/
-
     /*--- Only apply thickened flame correction factor to sources for steady problems. ---*/
-    su2double F_source = config->GetTime_Domain() ? F_default : 1.0 / F;
     unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector, F_source);
 
     if (ignition) {
@@ -117,8 +136,10 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
                 spark_radius = flamelet_config_options.spark_init[3];
       dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), flamelet_config_options.spark_init.data());
       if (dist_from_center < pow(spark_radius,2)) {
-        for (auto iVar = 0u; iVar < nVar; iVar++)
-          nodes->SetScalarSource(i_point, iVar,  nodes->GetScalarSources(i_point)[iVar] + flamelet_config_options.spark_reaction_rates[iVar] / F);
+        for (auto iVar = 0u; iVar < nVar; iVar++) {
+          su2double source_total = nodes->GetScalarSources(i_point)[iVar] + F_source * flamelet_config_options.spark_reaction_rates[iVar];
+          nodes->SetScalarSource(i_point, iVar,  source_total);
+        }
       }
     }
 
@@ -822,7 +843,50 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo
 
 
 su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) {
-  su2double flame_vol = pow(flamelet_config_options.flame_lengthscale,nDim);
-  su2double F = max(1.0, geometry->nodes->GetVolume(iPoint) / flame_vol);
+  su2double F{1.0};
+  if (global_flame_thickness != default_flame_thickness) {
+    su2double max_flame_vol = pow(global_flame_thickness/extra_flame_thickness_correction, nDim);
+    F = max(1.0, geometry->nodes->GetVolume(iPoint) / max_flame_vol);
+  }
   return F;
+}
+
+su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry * geometry, CSolver **solver_container) {
+  CFlowVariable* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
+  su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{0.0}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{0.0}, Tmax_global{0.0};
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
+      su2double pv_local = nodes->GetSolution(iPoint, I_PROGVAR);
+      su2double wdist = geometry->nodes->GetWall_Distance(iPoint);
+      if (global_flame_thickness != default_flame_thickness)  {
+        su2double T_local = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint);
+        Tmax_local = max(T_local, Tmax_local);
+        
+        pvmax_local = max(pv_local, pvmax_local);
+        pvmin_local = min(pv_local, pvmin_local);
+        su2double gradpv{0.0};
+        su2double grad_pv{0.0};
+        
+        su2double fac = min(wdist / global_flame_thickness, 1.0);
+        for (auto iDim=0u; iDim<nDim; iDim++) {
+          su2double dpvdx = nodes->GetGradient(iPoint, I_PROGVAR, iDim);
+          
+          gradpv+= pow(dpvdx,2);
+        };
+        gradpv_local = max(gradpv_local, fac*sqrt(gradpv));
+      
+      }
+    
+  }
+  END_SU2_OMP_FOR
+  SU2_MPI::Barrier(SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&gradpv_local, &gradpv_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&pvmax_local, &pvmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&pvmin_local, &pvmin_global, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&Tmax_local, &Tmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+  su2double flame_thickness{default_flame_thickness};
+  if (Tmax_global > 800.0) flame_thickness = (pvmax_global - pvmin_global) / gradpv_global;
+
+  
+  return flame_thickness;
 }
\ No newline at end of file

From d9e5b59dabf016667392d71abd2a9fd16b8fd1c0 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Thu, 19 Mar 2026 11:08:38 +0100
Subject: [PATCH 09/23] Added default flame thickness option

---
 Common/src/CConfig.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 33eeaec75f6..9e0ecf5a8ef 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1430,7 +1430,7 @@ void CConfig::SetConfig_Options() {
   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);
+  addDoubleOption("FLAME_LENGTHSCALE", flamelet_ParsedOptions.flame_lengthscale, flamelet_ParsedOptions.default_lengthscale);
 
   /*--- Options related to mass diffusivity and thereby the species solver. ---*/
 
@@ -5794,7 +5794,7 @@ 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;
 
-    if (flamelet_ParsedOptions.flame_lengthscale <= 0.0)
+    if (flamelet_ParsedOptions.flame_lengthscale <= 0.0 && flamelet_ParsedOptions.flame_lengthscale != flamelet_ParsedOptions.default_lengthscale)
       SU2_MPI::Error("Flame length scale value should be positive.", CURRENT_FUNCTION);
 
   }

From 6963939d310e6a854a0be0e3fae3eaadd5236024 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Fri, 20 Mar 2026 15:38:02 +0100
Subject: [PATCH 10/23] Removed flame lengthscale option and automatically
 approximate flame thickness

---
 Common/include/option_structure.hpp           |   3 +-
 Common/src/CConfig.cpp                        |   8 +-
 .../solvers/CSpeciesFlameletSolver.hpp        |  19 +++-
 .../src/solvers/CSpeciesFlameletSolver.cpp    | 102 ++++++++++--------
 config_template.cfg                           |   5 +
 5 files changed, 81 insertions(+), 56 deletions(-)

diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index dff478254dd..7263b98acc3 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -1446,8 +1446,7 @@ 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.*/
+  bool thickenedflame_correction{true}; /*!< \brief Thickened flame correction. */
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index f3d5a3bd956..de803a1640b 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1429,8 +1429,8 @@ 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, flamelet_ParsedOptions.default_lengthscale);
+  /*!\brief THICKENED_FLAME_CORRECTION \n DESCRIPTION: Dampen source terms and enhance diffusion based on the flame length scale. \ingroup Config*/
+  addBoolOption("THICKENED_FLAME_CORRECTION", flamelet_ParsedOptions.thickenedflame_correction, true);
 
   /*--- Options related to mass diffusivity and thereby the species solver. ---*/
 
@@ -5795,10 +5795,6 @@ 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 && flamelet_ParsedOptions.flame_lengthscale != flamelet_ParsedOptions.default_lengthscale)
-      SU2_MPI::Error("Flame length scale value should be positive.", CURRENT_FUNCTION);
-
   }
 
   if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) {
diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
index 26f44286691..1da87f4e878 100644
--- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
+++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
@@ -37,6 +37,9 @@
  */
 class CSpeciesFlameletSolver final : public CSpeciesSolver {
  private:
+  const su2double default_flame_thickness{1.0};
+  su2double global_flame_thickness;
+  bool calc_flame_thickness{false};
   FluidFlamelet_ParsedOptions flamelet_config_options;
   bool include_mixture_fraction = false; /*!< \brief include mixture fraction as a controlling variable. */
   /*!
@@ -113,7 +116,15 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \param[in] iPoint - node ID.
    * \return - flame thickness correction factor.
    */                                              
-  su2double ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint);   
+  su2double ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) const;   
+
+  /*!
+   * \brief Approximate the minimum flame thickness value used for the thickened flame model.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \return - approximate flame thickness value.
+   */
+  su2double GetOverallFlameThickness(CGeometry * geometry,  CSolver **solver_container) const;
 
  public:
   /*!
@@ -204,4 +215,10 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    */
   void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics,
                         const CConfig* config) final;
+
+  /*!
+   * \brief Obtain the overall flame thickness value.
+   * \return flame thickness value.
+   */
+  virtual su2double GetFlameThickness() const override {return global_flame_thickness;}
 };
diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index f5406aa5189..015eac3daf6 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -39,14 +39,8 @@ CSpeciesFlameletSolver::CSpeciesFlameletSolver(CGeometry* geometry, CConfig* con
 
   /*--- Retrieve options from config. ---*/
   flamelet_config_options = config->GetFlameletParsedOptions();
-
-  if (flamelet_config_options.flame_lengthscale != flamelet_config_options.default_lengthscale) {
-    global_flame_thickness = flamelet_config_options.flame_lengthscale;
-    calc_flame_thickness = false;
-  } else {
-    global_flame_thickness = default_flame_thickness;
-    calc_flame_thickness = true;
-  }
+  global_flame_thickness = default_flame_thickness;
+  calc_flame_thickness = flamelet_config_options.thickenedflame_correction;
 
   /*--- Dimension of the problem. ---*/
   nVar = flamelet_config_options.n_scalars;
@@ -74,12 +68,8 @@ CSpeciesFlameletSolver::CSpeciesFlameletSolver(CGeometry* geometry, CConfig* con
   /*--- Add the solver name. ---*/
   SolverName = "FLAMELET";
 
-  if (rank==MASTER_NODE) {
-    if (calc_flame_thickness) {
-      cout << "Lengthscale for thickened flame model is automatically determined." << endl;
-    } else {
-      cout << "Thickened flame model is applied in cells larger than " << global_flame_thickness*1000 << " mm" << endl;
-    }
+  if (calc_flame_thickness && rank==MASTER_NODE) {
+    cout << "Applying thickened flame source and diffusion correction." << endl;
   }
 }
 
@@ -111,19 +101,29 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
   }
   SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);)
 
-  su2double test_thickness = GetOverallFlameThickness(geometry, solver_container);
-  if (test_thickness < global_flame_thickness) {
-    global_flame_thickness = test_thickness;
-  } 
-  if (rank==MASTER_NODE) cout << global_flame_thickness << endl;
+  /* Update global flame thickness value. */
+  if (calc_flame_thickness) {
+    su2double test_thickness = GetOverallFlameThickness(geometry, solver_container);
+    if (test_thickness < global_flame_thickness) {
+      global_flame_thickness = test_thickness;
+    } else {
+      global_flame_thickness = 0.95*global_flame_thickness + 0.05*test_thickness;
+    }
+  }
+
+  /* Flame thickness correction factors */
+  su2double F{1.0}, F_source{1.0};
+
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto i_point = 0u; i_point < nPoint; i_point++) {
     CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel();
     su2double* scalars = nodes->GetSolution(i_point);
 
     /*--- Calculate correction factor for flame propagation on coarse grids. ---*/
-    const su2double F = ThickenedFlameCorrection(geometry, i_point);
-    const su2double F_source = 1.0 / F;
+    if (calc_flame_thickness) {
+      F = ThickenedFlameCorrection(geometry, i_point);
+      F_source = 1.0 / F;
+    }
 
     for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar];
 
@@ -842,51 +842,59 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo
 }
 
 
-su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) {
+su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) const {
   su2double F{1.0};
   if (global_flame_thickness != default_flame_thickness) {
-    su2double max_flame_vol = pow(global_flame_thickness/extra_flame_thickness_correction, nDim);
+    su2double max_flame_vol = pow(global_flame_thickness, nDim);
     F = max(1.0, geometry->nodes->GetVolume(iPoint) / max_flame_vol);
   }
   return F;
 }
 
-su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry * geometry, CSolver **solver_container) {
-  CFlowVariable* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
-  su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{0.0}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{0.0}, Tmax_global{0.0};
+su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry * geometry, CSolver **solver_container) const {
+  CFlowVariable const * flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
+  su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{0.0}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{-1e6},Tmax_global{0.0};
+  
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
       su2double pv_local = nodes->GetSolution(iPoint, I_PROGVAR);
-      su2double wdist = geometry->nodes->GetWall_Distance(iPoint);
-      if (global_flame_thickness != default_flame_thickness)  {
-        su2double T_local = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint);
-        Tmax_local = max(T_local, Tmax_local);
-        
-        pvmax_local = max(pv_local, pvmax_local);
-        pvmin_local = min(pv_local, pvmin_local);
-        su2double gradpv{0.0};
-        su2double grad_pv{0.0};
-        
-        su2double fac = min(wdist / global_flame_thickness, 1.0);
-        for (auto iDim=0u; iDim<nDim; iDim++) {
-          su2double dpvdx = nodes->GetGradient(iPoint, I_PROGVAR, iDim);
-          
-          gradpv+= pow(dpvdx,2);
-        };
-        gradpv_local = max(gradpv_local, fac*sqrt(gradpv));
+      su2double T_local = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint);
       
+      /* Parallel projection of progress variable gradient against velocity */
+      su2double proj_grad_pv_u[MAXNDIM];
+      for (auto iDim=0u; iDim < nDim; iDim++) {
+        su2double gradpv = nodes->GetGradient(iPoint, I_PROGVAR, iDim);
+        su2double val_u = flowNodes->GetVelocity(iPoint, iDim);
+        proj_grad_pv_u[iDim] = gradpv * val_u * val_u / (flowNodes->GetVelocity2(iPoint) + 1e-5);
       }
-    
+
+      /* Parallel projection of temperature gradient against projected progress variable gradient */
+      su2double proj_grad_T_u = 0, gradT2{0}, mag_gradT{0};
+      for (auto iDim=0u; iDim < nDim; iDim++) {
+        su2double gradT = flowNodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim);
+        proj_grad_T_u += gradT * proj_grad_pv_u[iDim];
+        gradT2 += gradT*gradT;
+      }
+      mag_gradT = sqrt(gradT2);
+      proj_grad_T_u /= (gradT2 +  1e-5);
+      proj_grad_T_u *= mag_gradT;
+
+      /* Update minimum and maximum values. */
+      gradpv_local = max(gradpv_local, proj_grad_T_u);
+      pvmax_local = max(pvmax_local, pv_local);
+      pvmin_local = min(pvmin_local, pv_local);
+      Tmax_local = max(Tmax_local, T_local);
   }
   END_SU2_OMP_FOR
   SU2_MPI::Barrier(SU2_MPI::GetComm());
   SU2_MPI::Allreduce(&gradpv_local, &gradpv_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
-  SU2_MPI::Allreduce(&pvmax_local, &pvmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
   SU2_MPI::Allreduce(&pvmin_local, &pvmin_global, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&pvmax_local, &pvmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
   SU2_MPI::Allreduce(&Tmax_local, &Tmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+
+  /* Update flame thickness value. */
   su2double flame_thickness{default_flame_thickness};
-  if (Tmax_global > 800.0) flame_thickness = (pvmax_global - pvmin_global) / gradpv_global;
+  if (Tmax_global > 800.0) flame_thickness = (pvmax_global - pvmin_global) / (gradpv_global+1e-5);
 
-  
   return flame_thickness;
 }
\ No newline at end of file
diff --git a/config_template.cfg b/config_template.cfg
index 0dd67764bdf..1201c5eb076 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -974,6 +974,11 @@ SPARK_INIT= (0.004, 0.0, 0.0, 1e-4, 100, 10)
 % The number of terms should equate the total number of species in the flamelet problem.
 SPARK_REACTION_RATES= (1000, 0, 0)
 
+% Thickened flame correction
+% Approximate the thickness of the flame and suppress source terms and enhance diffusivity in cells with a lengthscale higher 
+% than the flame thickness.
+THICKENED_FLAME_CORRECTION= YES
+
 %
 % --------------------- INVERSE DESIGN SIMULATION -----------------------------%
 %

From e8b31f95aa3b2d50d0e92ba914a429429de712bb Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Fri, 20 Mar 2026 15:38:25 +0100
Subject: [PATCH 11/23] Include flame thickness in convergence history

---
 SU2_CFD/include/solvers/CSolver.hpp | 6 ++++++
 SU2_CFD/src/output/CFlowOutput.cpp  | 5 +++++
 2 files changed, 11 insertions(+)

diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
index 67a8bb8112c..da5244c9298 100644
--- a/SU2_CFD/include/solvers/CSolver.hpp
+++ b/SU2_CFD/include/solvers/CSolver.hpp
@@ -580,6 +580,12 @@ class CSolver {
    */
   inline virtual void SetPrimitive_Limiter(CGeometry *geometry, const CConfig *config) { }
 
+  /*!
+   * \brief A virtual member.
+   * \return flame thickness value.
+   */
+  virtual su2double GetFlameThickness() const {return 1.0;}
+
   /*!
    * \brief Compute the projection of a variable for MUSCL reconstruction.
    * \note The result should be halved when added to i (or subtracted from j).
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 446db73dc13..2ca35e74c7c 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1069,6 +1069,8 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) {
         const auto& CV_name = flamelet_config_options.controlling_variable_names[iCV];
         AddHistoryOutput("RMS_"+CV_name, "rms["+CV_name+"]",ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean squared residual of " + CV_name + " controlling variable equation.", HistoryFieldType::RESIDUAL);
       }
+      if (flamelet_config_options.thickenedflame_correction)
+        AddHistoryOutput("THICKNESS","flamethickness",ScreenOutputFormat::FIXED, "RMS_RES", "Flame thickness used for thickened flame correction model.", HistoryFieldType::RESIDUAL);
 
       /*--- auxiliary species transport ---*/
       for (auto i_scalar = 0u; i_scalar < flamelet_config_options.n_user_scalars; i_scalar++){
@@ -1304,6 +1306,9 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co
         }
       }
 
+      if (flamelet_config_options.thickenedflame_correction)
+        SetHistoryOutputValue("THICKNESS", solver[SPECIES_SOL]->GetFlameThickness());
+
       SetHistoryOutputValue("LINSOL_ITER_FLAMELET", solver[SPECIES_SOL]->GetIterLinSolver());
       SetHistoryOutputValue("LINSOL_RESIDUAL_FLAMELET", log10(solver[SPECIES_SOL]->GetResLinSolver()));
     }

From 5dee84805c872b929869c291b08c66b489e1437f Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Fri, 20 Mar 2026 15:39:17 +0100
Subject: [PATCH 12/23] precommit formatting

---
 config_template.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/config_template.cfg b/config_template.cfg
index 1201c5eb076..62ca6b4e464 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -975,7 +975,7 @@ SPARK_INIT= (0.004, 0.0, 0.0, 1e-4, 100, 10)
 SPARK_REACTION_RATES= (1000, 0, 0)
 
 % Thickened flame correction
-% Approximate the thickness of the flame and suppress source terms and enhance diffusivity in cells with a lengthscale higher 
+% Approximate the thickness of the flame and suppress source terms and enhance diffusivity in cells with a lengthscale higher
 % than the flame thickness.
 THICKENED_FLAME_CORRECTION= YES
 

From 101042975c202253fa91a6f66cb68c3fe253b438 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Fri, 20 Mar 2026 15:40:21 +0100
Subject: [PATCH 13/23] code factor

---
 SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
index 1da87f4e878..a882068cd59 100644
--- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
+++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
@@ -220,5 +220,5 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \brief Obtain the overall flame thickness value.
    * \return flame thickness value.
    */
-  virtual su2double GetFlameThickness() const override {return global_flame_thickness;}
+  virtual su2double GetFlameThickness() const {return global_flame_thickness;}
 };

From 98de6ac28856914a14db11fec84937eff59dd4cb Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 09:33:42 +0100
Subject: [PATCH 14/23] Applied GitHub code suggestion

---
 SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 015eac3daf6..23bfe5534b4 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -844,7 +844,7 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo
 
 su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) const {
   su2double F{1.0};
-  if (global_flame_thickness != default_flame_thickness) {
+  if (fabs(global_flame_thickness - default_flame_thickness) > EPS * max(1.0, fabs(default_flame_thickness))) {
     su2double max_flame_vol = pow(global_flame_thickness, nDim);
     F = max(1.0, geometry->nodes->GetVolume(iPoint) / max_flame_vol);
   }

From 4ff73c6556ea015eb2c1fb8ce9f42cad1bcb3c02 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 09:36:04 +0100
Subject: [PATCH 15/23] Improved description of  coarse grid correction

---
 Common/src/CConfig.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index de803a1640b..f780df917db 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1429,7 +1429,7 @@ 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 THICKENED_FLAME_CORRECTION \n DESCRIPTION: Dampen source terms and enhance diffusion based on the flame length scale. \ingroup Config*/
+  /*!\brief THICKENED_FLAME_CORRECTION \n DESCRIPTION: Coarse grid correction for source terms and diffusive fluxes in reacting flows. \ingroup Config*/
   addBoolOption("THICKENED_FLAME_CORRECTION", flamelet_ParsedOptions.thickenedflame_correction, true);
 
   /*--- Options related to mass diffusivity and thereby the species solver. ---*/

From 746a543ede4b0729bb3778eb69e1bbc0515c2b5f Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 09:41:10 +0100
Subject: [PATCH 16/23] Formatting suggestions Pedro

---
 SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp | 6 +++---
 SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp     | 5 ++---
 2 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
index a882068cd59..a792bba6e80 100644
--- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
+++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp
@@ -116,7 +116,7 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \param[in] iPoint - node ID.
    * \return - flame thickness correction factor.
    */                                              
-  su2double ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) const;   
+  su2double ThickenedFlameCorrection(const CGeometry* geometry, unsigned long iPoint) const;   
 
   /*!
    * \brief Approximate the minimum flame thickness value used for the thickened flame model.
@@ -124,7 +124,7 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \param[in] solver_container - Container vector with all the solutions.
    * \return - approximate flame thickness value.
    */
-  su2double GetOverallFlameThickness(CGeometry * geometry,  CSolver **solver_container) const;
+  su2double GetOverallFlameThickness(CGeometry* geometry,  CSolver** solver_container) const;
 
  public:
   /*!
@@ -220,5 +220,5 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \brief Obtain the overall flame thickness value.
    * \return flame thickness value.
    */
-  virtual su2double GetFlameThickness() const {return global_flame_thickness;}
+  su2double GetFlameThickness() const override {return global_flame_thickness;}
 };
diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 23bfe5534b4..537ab80eb67 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -418,7 +418,6 @@ void CSpeciesFlameletSolver::SetPreconditioner(CGeometry* geometry, CSolver** so
 
 void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solver_container,
                                              CNumerics** numerics_container, CConfig* config, unsigned short iMesh) {
-   
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto i_point = 0u; i_point < nPointDomain; i_point++) {
     /*--- Add source terms from the lookup table directly to the residual. ---*/
@@ -842,7 +841,7 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo
 }
 
 
-su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) const {
+su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(const CGeometry* geometry, unsigned long iPoint) const {
   su2double F{1.0};
   if (fabs(global_flame_thickness - default_flame_thickness) > EPS * max(1.0, fabs(default_flame_thickness))) {
     su2double max_flame_vol = pow(global_flame_thickness, nDim);
@@ -851,7 +850,7 @@ su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geo
   return F;
 }
 
-su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry * geometry, CSolver **solver_container) const {
+su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry, CSolver** solver_container) const {
   CFlowVariable const * flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
   su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{0.0}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{-1e6},Tmax_global{0.0};
   

From 7cdd149264821a546221a8071186e5ea9b558c44 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 09:49:35 +0100
Subject: [PATCH 17/23] Flame thickness in history output now defined as a
 user-defined coefficient under its own group

---
 SU2_CFD/src/output/CFlowOutput.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 2ca35e74c7c..5c0d44ef5f6 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1070,7 +1070,7 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) {
         AddHistoryOutput("RMS_"+CV_name, "rms["+CV_name+"]",ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean squared residual of " + CV_name + " controlling variable equation.", HistoryFieldType::RESIDUAL);
       }
       if (flamelet_config_options.thickenedflame_correction)
-        AddHistoryOutput("THICKNESS","flamethickness",ScreenOutputFormat::FIXED, "RMS_RES", "Flame thickness used for thickened flame correction model.", HistoryFieldType::RESIDUAL);
+        AddHistoryOutput("THICKNESS","flamethickness",ScreenOutputFormat::FIXED, "FLAME_THICKNESS", "Flame thickness used for thickened flame correction model.", HistoryFieldType::COEFFICIENT);
 
       /*--- auxiliary species transport ---*/
       for (auto i_scalar = 0u; i_scalar < flamelet_config_options.n_user_scalars; i_scalar++){

From 7015adcaafca35a94ccf750b694ba14b6e165f23 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 09:52:26 +0100
Subject: [PATCH 18/23] const correction

---
 SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 537ab80eb67..0d476747e17 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -851,7 +851,7 @@ su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(const CGeometry* geom
 }
 
 su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry, CSolver** solver_container) const {
-  CFlowVariable const * flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
+  const CFlowVariable* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
   su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{0.0}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{-1e6},Tmax_global{0.0};
   
   SU2_OMP_FOR_STAT(omp_chunk_size)

From 671c8c384f1b362d49a917529a680b5855b33ce2 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 10:07:39 +0100
Subject: [PATCH 19/23] Single allreduce of four variables instead of four
 separate allreduce calls.

---
 .../src/solvers/CSpeciesFlameletSolver.cpp    | 31 +++++++++++++------
 1 file changed, 21 insertions(+), 10 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 0d476747e17..9080bfbc450 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -852,10 +852,10 @@ su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(const CGeometry* geom
 
 su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry, CSolver** solver_container) const {
   const CFlowVariable* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
-  su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{0.0}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{-1e6},Tmax_global{0.0};
+  su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{-1e3}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{-1e6},Tmax_global{0.0};
   
   SU2_OMP_FOR_STAT(omp_chunk_size)
-  for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
+  for (auto iPoint = 0u; iPoint < nPointDomain; iPoint++) {
       su2double pv_local = nodes->GetSolution(iPoint, I_PROGVAR);
       su2double T_local = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint);
       
@@ -875,25 +875,36 @@ su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry,
         gradT2 += gradT*gradT;
       }
       mag_gradT = sqrt(gradT2);
-      proj_grad_T_u /= (gradT2 +  1e-5);
+      proj_grad_T_u /= (gradT2 +  EPS);
       proj_grad_T_u *= mag_gradT;
 
       /* Update minimum and maximum values. */
       gradpv_local = max(gradpv_local, proj_grad_T_u);
       pvmax_local = max(pvmax_local, pv_local);
-      pvmin_local = min(pvmin_local, pv_local);
+      pvmin_local = max(pvmin_local, -pv_local);
       Tmax_local = max(Tmax_local, T_local);
   }
   END_SU2_OMP_FOR
   SU2_MPI::Barrier(SU2_MPI::GetComm());
-  SU2_MPI::Allreduce(&gradpv_local, &gradpv_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
-  SU2_MPI::Allreduce(&pvmin_local, &pvmin_global, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm());
-  SU2_MPI::Allreduce(&pvmax_local, &pvmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
-  SU2_MPI::Allreduce(&Tmax_local, &Tmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
-
+  su2double* MyFlameThickness = new su2double[4];
+  su2double* TotalFlameThickness = new su2double[4]; 
+  MyFlameThickness[0] = gradpv_local;
+  MyFlameThickness[1] = pvmin_local;
+  MyFlameThickness[2] = pvmax_local;
+  MyFlameThickness[3] = Tmax_local;
+  
+  SU2_MPI::Allreduce(MyFlameThickness, TotalFlameThickness, 4, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+  gradpv_global = TotalFlameThickness[0];
+  pvmin_global = -TotalFlameThickness[1];
+  pvmax_global = TotalFlameThickness[2];
+  Tmax_global = TotalFlameThickness[3];
+
+  delete [] MyFlameThickness;
+  delete [] TotalFlameThickness;
+  
   /* Update flame thickness value. */
   su2double flame_thickness{default_flame_thickness};
-  if (Tmax_global > 800.0) flame_thickness = (pvmax_global - pvmin_global) / (gradpv_global+1e-5);
+  if (Tmax_global > 800.0) flame_thickness = (pvmax_global - pvmin_global) / (gradpv_global+EPS);
 
   return flame_thickness;
 }
\ No newline at end of file

From 566618e71a17b104109c7c94638b949d4f809a77 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 10:21:47 +0100
Subject: [PATCH 20/23] Geometrytoolbox for gradient projection

---
 .../src/solvers/CSpeciesFlameletSolver.cpp    | 21 +++++++++----------
 1 file changed, 10 insertions(+), 11 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 9080bfbc450..8f471f7bb47 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -860,22 +860,22 @@ su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry,
       su2double T_local = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint);
       
       /* Parallel projection of progress variable gradient against velocity */
-      su2double proj_grad_pv_u[MAXNDIM];
+      su2double proj_grad_pv_u[MAXNDIM]={0};
       for (auto iDim=0u; iDim < nDim; iDim++) {
         su2double gradpv = nodes->GetGradient(iPoint, I_PROGVAR, iDim);
         su2double val_u = flowNodes->GetVelocity(iPoint, iDim);
-        proj_grad_pv_u[iDim] = gradpv * val_u * val_u / (flowNodes->GetVelocity2(iPoint) + 1e-5);
+        proj_grad_pv_u[iDim] = gradpv * val_u * val_u / (flowNodes->GetVelocity2(iPoint) + EPS);
       }
 
       /* Parallel projection of temperature gradient against projected progress variable gradient */
-      su2double proj_grad_T_u = 0, gradT2{0}, mag_gradT{0};
-      for (auto iDim=0u; iDim < nDim; iDim++) {
-        su2double gradT = flowNodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim);
-        proj_grad_T_u += gradT * proj_grad_pv_u[iDim];
-        gradT2 += gradT*gradT;
-      }
-      mag_gradT = sqrt(gradT2);
-      proj_grad_T_u /= (gradT2 +  EPS);
+      su2double gradT[MAXNDIM]={0};
+      for (auto iDim=0u; iDim < nDim; iDim++) 
+        gradT[iDim] = flowNodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim);
+      
+      su2double proj_grad_T_u = GeometryToolbox::DotProduct(nDim, gradT, proj_grad_pv_u);
+      su2double mag_gradT = GeometryToolbox::Norm(nDim, gradT);
+      
+      proj_grad_T_u /= (pow(mag_gradT,2) +  EPS);
       proj_grad_T_u *= mag_gradT;
 
       /* Update minimum and maximum values. */
@@ -885,7 +885,6 @@ su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry,
       Tmax_local = max(Tmax_local, T_local);
   }
   END_SU2_OMP_FOR
-  SU2_MPI::Barrier(SU2_MPI::GetComm());
   su2double* MyFlameThickness = new su2double[4];
   su2double* TotalFlameThickness = new su2double[4]; 
   MyFlameThickness[0] = gradpv_local;

From 8f1f2a170f20861d8da23254662f14bedc1b7318 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 11:24:16 +0100
Subject: [PATCH 21/23] Added ignition temperature option

---
 Common/include/option_structure.hpp | 1 +
 Common/src/CConfig.cpp              | 3 +++
 config_template.cfg                 | 6 +++++-
 3 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 7263b98acc3..b524a8b583b 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -1447,6 +1447,7 @@ struct FluidFlamelet_ParsedOptions {
   unsigned short nspark;           /*!< \brief Number of source terms for spark initialization. */
   bool preferential_diffusion = false;  /*!< \brief Preferential diffusion physics for flamelet solver.*/
   bool thickenedflame_correction{true}; /*!< \brief Thickened flame correction. */
+  su2double ignition_temperature = 1100.0; /*!< \brief Temperature above which the solution is considered ignited. */
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index f780df917db..911bd6e3b20 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1432,6 +1432,9 @@ void CConfig::SetConfig_Options() {
   /*!\brief THICKENED_FLAME_CORRECTION \n DESCRIPTION: Coarse grid correction for source terms and diffusive fluxes in reacting flows. \ingroup Config*/
   addBoolOption("THICKENED_FLAME_CORRECTION", flamelet_ParsedOptions.thickenedflame_correction, true);
 
+  /*!\brief IGNITION_TEMPERATURE \n DESCRIPTION: Temperature above which the solution is considered ignited. */
+  addDoubleOption("IGNITION_TEMPERATURE", flamelet_ParsedOptions.ignition_temperature, 1100.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*/
diff --git a/config_template.cfg b/config_template.cfg
index 62ca6b4e464..1c6ff7b1f9c 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -934,7 +934,7 @@ CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL)
 % Method used to ignite the solution:
 % FLAME_FRONT : the flame is initialized using a plane, defined by a point and a normal. On one side, the solution is initialized
 % using 'burnt' conditions and on the other side 'unburnt' conditions. The normal points in the direction of the 'burnt'
-% condition.
+% condition. The 'burnt' condition is defined as the value of the progress variable for which the temperature exceeds the IGNITION_TEMPERATURE.
 % SPARK : the solution is ignited through application of a set of source terms within a specified region for a set number
 % of solver iterations. This artificial spark can be applied after a certain number of iterations has passed, allowing for
 % the flow to evolve before igniting the mixture.
@@ -979,6 +979,10 @@ SPARK_REACTION_RATES= (1000, 0, 0)
 % than the flame thickness.
 THICKENED_FLAME_CORRECTION= YES
 
+% Flame ignition temperature
+% Flame thickness calculations are performed when the maximum temperature in the flow solution exceeds this value. 
+IGNITION_TEMPERATURE= 1100.0
+
 %
 % --------------------- INVERSE DESIGN SIMULATION -----------------------------%
 %

From d6ddccb6bc393904f50e1b5a3cc30c3f25ca2a02 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 11:24:43 +0100
Subject: [PATCH 22/23] Clip global flame thickness to one meter maximum.

---
 .../src/solvers/CSpeciesFlameletSolver.cpp    | 28 +++++++++----------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
index 8f471f7bb47..6e121e6710c 100644
--- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
@@ -103,7 +103,8 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
 
   /* Update global flame thickness value. */
   if (calc_flame_thickness) {
-    su2double test_thickness = GetOverallFlameThickness(geometry, solver_container);
+    su2double calc_thickness = GetOverallFlameThickness(geometry, solver_container);
+    su2double test_thickness = min(default_flame_thickness, calc_thickness);
     if (test_thickness < global_flame_thickness) {
       global_flame_thickness = test_thickness;
     } else {
@@ -833,7 +834,7 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo
   bool outside = false;
   while (!outside) {
     fluid_model->SetTDState_T(300, scalars);
-    if (fluid_model->GetExtrapolation() == 1 || (fluid_model->GetTemperature()>1000.)) outside = true;
+    if (fluid_model->GetExtrapolation() == 1 || fluid_model->GetTemperature()>flamelet_config_options.ignition_temperature) outside = true;
     scalars[I_PROGVAR] += delta;
   }
   su2double pv_burnt = scalars[I_PROGVAR] - delta;
@@ -852,7 +853,7 @@ su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(const CGeometry* geom
 
 su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry, CSolver** solver_container) const {
   const CFlowVariable* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
-  su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{-1e3}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{-1e6},Tmax_global{0.0};
+  su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{1e3}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{-1e6},Tmax_global{0.0};
   
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto iPoint = 0u; iPoint < nPointDomain; iPoint++) {
@@ -881,29 +882,28 @@ su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry,
       /* Update minimum and maximum values. */
       gradpv_local = max(gradpv_local, proj_grad_T_u);
       pvmax_local = max(pvmax_local, pv_local);
-      pvmin_local = max(pvmin_local, -pv_local);
+      pvmin_local = min(pvmin_local, pv_local);
       Tmax_local = max(Tmax_local, T_local);
   }
   END_SU2_OMP_FOR
-  su2double* MyFlameThickness = new su2double[4];
-  su2double* TotalFlameThickness = new su2double[4]; 
+  su2double* MyFlameThickness = new su2double[3];
+  su2double* TotalFlameThickness = new su2double[3]; 
   MyFlameThickness[0] = gradpv_local;
-  MyFlameThickness[1] = pvmin_local;
-  MyFlameThickness[2] = pvmax_local;
-  MyFlameThickness[3] = Tmax_local;
+  MyFlameThickness[1] = pvmax_local;
+  MyFlameThickness[2] = Tmax_local;
   
-  SU2_MPI::Allreduce(MyFlameThickness, TotalFlameThickness, 4, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(MyFlameThickness, TotalFlameThickness, 3, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&pvmin_local, &pvmin_global, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm());
   gradpv_global = TotalFlameThickness[0];
-  pvmin_global = -TotalFlameThickness[1];
-  pvmax_global = TotalFlameThickness[2];
-  Tmax_global = TotalFlameThickness[3];
+  pvmax_global = TotalFlameThickness[1];
+  Tmax_global = TotalFlameThickness[2];
 
   delete [] MyFlameThickness;
   delete [] TotalFlameThickness;
   
   /* Update flame thickness value. */
   su2double flame_thickness{default_flame_thickness};
-  if (Tmax_global > 800.0) flame_thickness = (pvmax_global - pvmin_global) / (gradpv_global+EPS);
+  if (Tmax_global > flamelet_config_options.ignition_temperature) flame_thickness = (pvmax_global - pvmin_global) / (gradpv_global+EPS);
 
   return flame_thickness;
 }
\ No newline at end of file

From e9d88dcc988feca261e825363448c794337a6575 Mon Sep 17 00:00:00 2001
From: EvertBunschoten <evertbun@xs4all.nl>
Date: Wed, 25 Mar 2026 11:27:45 +0100
Subject: [PATCH 23/23] precommit

---
 config_template.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/config_template.cfg b/config_template.cfg
index 1c6ff7b1f9c..0ecf0d5e82d 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -980,7 +980,7 @@ SPARK_REACTION_RATES= (1000, 0, 0)
 THICKENED_FLAME_CORRECTION= YES
 
 % Flame ignition temperature
-% Flame thickness calculations are performed when the maximum temperature in the flow solution exceeds this value. 
+% Flame thickness calculations are performed when the maximum temperature in the flow solution exceeds this value.
 IGNITION_TEMPERATURE= 1100.0
 
 %