From 11a88fb506dcd8561817a7ea79cd461084eda637 Mon Sep 17 00:00:00 2001 From: EvertBunschoten Date: Tue, 17 Feb 2026 14:12:49 +0100 Subject: [PATCH 1/6] 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 Date: Tue, 17 Feb 2026 14:27:04 +0100 Subject: [PATCH 2/6] 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 Date: Fri, 20 Feb 2026 10:53:40 +0100 Subject: [PATCH 3/6] 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& scalars); + const vector& 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& 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 Date: Fri, 20 Feb 2026 10:55:02 +0100 Subject: [PATCH 4/6] 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& scalars) { + unsigned long iPoint, const vector& scalars, const su2double F) { /*--- Compute total source terms from the production and consumption. ---*/ vector 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 Date: Fri, 20 Feb 2026 13:17:38 +0100 Subject: [PATCH 5/6] 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 Date: Fri, 20 Feb 2026 13:43:49 +0100 Subject: [PATCH 6/6] 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 scalars_vector(nVar); + const su2double F_default{1.0}; unsigned long spark_iter_start, spark_duration; bool ignition = false; auto* flowNodes = su2staticcast_p(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; }