Fix flame initialization

.github/workflows/regression.yml

@@ -211,7 +211,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c fix_flame_initialization -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
@@ -260,7 +260,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests"
+          args: -b ${{github.ref}} -t develop -c fix_flame_initialization -s ${{matrix.testscript}} -a "--tapetests"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
@@ -308,7 +308,7 @@ jobs:
           PMIX_MCA_gds: hash
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan"
+          args: -b ${{github.ref}} -t develop -c fix_flame_initialization -s ${{matrix.testscript}} -a "--tsan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:260405-0054
         with:
@@ -353,7 +353,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan"
+          args: -b ${{github.ref}} -t develop -c fix_flame_initialization -s ${{matrix.testscript}} -a "--asan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054
         with:

Common/include/option_structure.hpp

@@ -1465,6 +1465,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_T_ignition = 5000;    /*!< \brief Ignition temperature for the flame, used for initialization. */
+
 };
 
 /*!

Common/src/CConfig.cpp

@@ -1419,6 +1419,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_INIT_IGNITION \n DESCRIPTION: Ignition temperature for the flame initialization \ingroup Config*/
+  addDoubleOption("FLAME_INIT_IGNITION", flamelet_ParsedOptions.Flame_T_ignition, 5000.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*/
@@ -5700,6 +5703,26 @@ 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;
+
+    /*--- Check that spark ignition has required parameters defined ---*/
+    if (flamelet_ParsedOptions.ignition_method == FLAMELET_INIT_TYPE::SPARK) {
+      /*--- Check if SPARK_INIT was explicitly set in config file ---*/
+      if (all_options.find("SPARK_INIT") != all_options.end()) {
+        SU2_MPI::Error("FLAME_INIT_METHOD= SPARK requires SPARK_INIT to be defined in the config file.", CURRENT_FUNCTION);
+      }
+      /*--- Check if SPARK_REACTION_RATES was explicitly set in config file ---*/
+      if (all_options.find("SPARK_REACTION_RATES") != all_options.end()) {
+        SU2_MPI::Error("FLAME_INIT_METHOD= SPARK requires SPARK_REACTION_RATES to be defined in the config file.", CURRENT_FUNCTION);
+      }
+      if (flamelet_ParsedOptions.nspark < flamelet_ParsedOptions.n_scalars) {
+        SU2_MPI::Error("SPARK_REACTION_RATES must have at least " + to_string(flamelet_ParsedOptions.n_scalars) +
+                       " values (one for each scalar variable), but only " + to_string(flamelet_ParsedOptions.nspark) + " were provided.", CURRENT_FUNCTION);
+      }
+    }
+    /*--- Check if flame ignition temperature is valid ---*/
+    if (flamelet_ParsedOptions.Flame_T_ignition <= Inc_Temperature_Init) {
+      SU2_MPI::Error("Flame ignition temperature must be higher than the initial temperature of the flow field.", CURRENT_FUNCTION);
+    }
   }
 
   if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) {

SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp

@@ -70,9 +70,10 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \brief Find maximum progress variable value within the manifold for the current solution.
    * \param[in] fluid_model - pointer to flamelet fluid model.
    * \param[in] scalars - local scalar solution.
+   * \param[in] T_ignition - ignition temperature - at this temperature, the progress variable is considered burnt.
    * \return - maximum progress variable value within manifold bounds.
    */
-  su2double GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalars);
+  su2double GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalars, const su2double T_ignition);
 
   /*!
    * \brief Retrieve scalar source terms from manifold.

SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp

@@ -103,9 +103,13 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
       su2double dist_from_center = 0,
                 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]);
+      su2double T_local = flowNodes->GetTemperature(i_point);
+      if (dist_from_center < pow(spark_radius,2) && T_local < flamelet_config_options.Flame_T_ignition) {
+        /*--- Add spark reaction rates to the sources that were just set by SetScalarSources ---*/
+        const su2double* current_sources = nodes->GetScalarSources(i_point);
+        for (auto iVar = 0u; iVar < nVar; iVar++) {
+          nodes->SetScalarSource(i_point, iVar, current_sources[iVar] + flamelet_config_options.spark_reaction_rates[iVar]);
+        }
       }
     }
 
@@ -184,7 +188,7 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver**
     su2double enth_inlet = config->GetSpecies_Init()[I_ENTH];
 
     su2double prog_burnt = 0, prog_unburnt, point_loc;
-    su2double scalar_init[MAXNVAR];
+    su2double scalar_init[MAXNVAR]= {0.0};
 
     if (rank == MASTER_NODE) {
       cout << "initial condition: T = " << temp_inlet << endl;
@@ -197,10 +201,11 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver**
           cout << "Ignition with a straight flame front" << endl;
           break;
         case FLAMELET_INIT_TYPE::SPARK:
-          cout << "Ignition with an artificial spark" << endl;
+          cout << "Ignition with an artificial spark at iteration "<< flamelet_config_options.spark_init[4]
+               << " for a duration of " << flamelet_config_options.spark_init[5] << " iterations." << endl;
           break;
         case FLAMELET_INIT_TYPE::NONE:
-          cout << "No solution ignition (cold flow)" << endl;
+          cout << "No solution ignition (cold flow or restart)" << endl;
           break;
         default:
           break;
@@ -216,14 +221,15 @@ 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();
-      if (flame_front_ignition) 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. ---*/
       n_not_iterated_local += GetEnthFromTemp(fluid_model_local, temp_inlet, config->GetSpecies_Init(), &enth_inlet);
       scalar_init[I_ENTH] = enth_inlet;
 
+      if (flame_front_ignition) prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init, flamelet_config_options.Flame_T_ignition);
+
       prog_unburnt = config->GetSpecies_Init()[I_PROGVAR];
       SU2_OMP_FOR_STAT(omp_chunk_size)
       for (unsigned long i_point = 0; i_point < nPoint; i_point++) {
@@ -805,16 +811,22 @@ unsigned long CSpeciesFlameletSolver::GetEnthFromTemp(CFluidModel* fluid_model,
   return exit_code;
 }
 
-su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalar_solution) {
+su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalar_solution, const su2double T_ignition) {
   SU2_ZONE_SCOPED
   su2double scalars[MAXNVAR], delta = 1e-3;
   for (auto iVar = 0u; iVar < nVar; iVar++) scalars[iVar] = scalar_solution[iVar];
   bool outside = false;
+  scalars[I_PROGVAR] += delta;
   while (!outside) {
-    fluid_model->SetTDState_T(300, scalars);
-    if (fluid_model->GetExtrapolation() == 1 || (fluid_model->GetTemperature()>1000.)) outside = true;
+    /*--- Note that 300.0 is a dummy temperature here and not used. ---*/
+    fluid_model->SetTDState_T(300.0, scalars);
+    if ((fluid_model->GetExtrapolation() == 1) || fluid_model->GetTemperature() > T_ignition) outside = true;
     scalars[I_PROGVAR] += delta;
   }
   su2double pv_burnt = scalars[I_PROGVAR] - delta;
+  if (rank == MASTER_NODE) {
+    cout << "Burnt progress variable determined from flamelet table: " << pv_burnt << endl;
+    cout << "Burnt temperature from flamelet table: " << fluid_model->GetTemperature() << endl;
+  }
   return pv_burnt;
 }

TestCases/flamelet/08_flame_init_method/flamelet_initialization.py

@@ -0,0 +1,144 @@
+#!/usr/bin/env python
+
+## \file run.py
+#  \brief Flame simulation initialization test case.
+#  \version 8.4.0 "Harrier"
+#
+# SU2 Project Website: https://su2code.github.io
+#
+# The SU2 Project is maintained by the SU2 Foundation
+# (http://su2foundation.org)
+#
+# Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
+#
+# SU2 is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# SU2 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+
+import pysu2 
+from mpi4py import MPI 
+
+# Generatl config options
+config_options = """
+FILENAMES_INTERPOLATOR= (MLP_Group1.mlp,MLP_Group2.mlp,MLP_Group3.mlp,MLP_Group4.mlp,MLP_Group5.mlp,MLP_Group6.mlp)
+
+INC_VELOCITY_INIT= (1e-5, 0.0, 0.0)
+MARKER_SYM=(x_minus, x_plus, y_plus, y_minus)
+
+SOLVER = INC_NAVIER_STOKES
+INC_DENSITY_MODEL= VARIABLE
+INC_ENERGY_EQUATION = YES
+INC_TEMPERATURE_INIT= 300
+INC_NONDIM= DIMENSIONAL
+
+
+FLUID_MODEL= FLUID_FLAMELET
+INTERPOLATION_METHOD= MLP
+
+CONTROLLING_VARIABLE_NAMES=(ProgressVariable,EnthalpyTot,MixtureFraction)
+CONTROLLING_VARIABLE_SOURCE_NAMES=(ProdRateTot_PV,NULL,NULL)
+
+KIND_SCALAR_MODEL= FLAMELET      
+DIFFUSIVITY_MODEL= FLAMELET
+VISCOSITY_MODEL= FLAMELET
+CONDUCTIVITY_MODEL= FLAMELET
+
+PREFERENTIAL_DIFFUSION=YES
+
+FLAME_INIT_METHOD= __FLAME_INIT_METHOD__
+SPARK_INIT=(5e-4, 2.5e-4, 0.0, 1e-4, 1, 50)
+SPARK_REACTION_RATES=(1.5e4, 0,0)
+FLAME_INIT= (4e-4, 0.00, 0.00, 1.0, 0.0, 0.0, 5e-5, 1.0)
+FLAME_INIT_IGNITION= __IGNITION_TEMP__
+
+CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR
+
+MUSCL_SPECIES= YES
+TIME_DISCRE_SPECIES= EULER_IMPLICIT
+SPECIES_CLIPPING= NO
+SPECIES_INIT=(-2.251135e-01,+2.262913e+03,+1.446752e-02)
+
+CFL_REDUCTION_SPECIES= 1.0
+
+INC_OUTLET_TYPE= PRESSURE_OUTLET
+
+NUM_METHOD_GRAD= GREEN_GAUSS
+CFL_NUMBER= 50
+CFL_ADAPT= NO
+ITER= __N_ITER__
+
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= LU_SGS
+LINEAR_SOLVER_ERROR= 1E-03
+LINEAR_SOLVER_ITER= 15
+
+CONV_NUM_METHOD_FLOW= FDS
+MUSCL_FLOW= YES
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+SCREEN_OUTPUT = CUSTOM
+HISTORY_OUTPUT = CUSTOM
+VOLUME_OUTPUT = SOLUTION
+
+MESH_FORMAT= RECTANGLE
+MESH_BOX_SIZE= ( 200, 50, 0 )
+MESH_BOX_LENGTH= ( 1e-3, 5e-4, 0 )
+
+OUTPUT_FILES=(RESTART)
+CUSTOM_OUTPUTS= 'Tprobe : Probe{TEMPERATURE}[5e-4, 2.5e-4];'
+
+"""
+
+comm = MPI.COMM_WORLD
+
+def main():
+
+    comm = MPI.COMM_WORLD
+
+    T_ignition = 1600
+    if comm.Get_rank() == 0:
+
+        # Prepare configuration files for flame front and spark ignition methods.
+        config_flame_front = config_options.replace("__FLAME_INIT_METHOD__", "FLAME_FRONT")
+        config_flame_front = config_flame_front.replace("__N_ITER__", "1")
+        config_flame_front = config_flame_front.replace("__IGNITION_TEMP__", "%f" % T_ignition)
+
+        config_spark = config_options.replace("__FLAME_INIT_METHOD__", "SPARK")
+        config_spark = config_spark.replace("__N_ITER__", "20")
+        config_spark = config_spark.replace("__IGNITION_TEMP__", "%f" % T_ignition)
+
+        with open("config_flame_front.cfg", "w") as f:
+            f.write(config_flame_front)
+        with open("config_spark.cfg", "w") as f:
+            f.write(config_spark)
+    comm.Barrier()
+
+    # Initialize simulation with flame front ignition
+    driver_flame_front = pysu2.CSinglezoneDriver("config_flame_front.cfg", 1, comm)
+    driver_flame_front.StartSolver()
+
+    # Retrieve temperature from the burnt zone of the initial flame front.
+    T_flame_front = driver_flame_front.GetOutputValue("Tprobe")
+    driver_flame_front.Finalize()
+
+    # Initialize simulation with spark ignition.
+    driver_spark = pysu2.CSinglezoneDriver("config_spark.cfg", 1, comm)
+    driver_spark.StartSolver()
+
+    # Retrieve temperature from center of the spark.
+    T_spark = driver_spark.GetOutputValue("Tprobe")
+    driver_spark.Finalize()
+    if comm.Get_rank() == 0:
+        print("| %i | %.2f | %.2f|" % (1, min(T_ignition, T_flame_front), min(T_ignition, T_spark)))
+
+if __name__ == '__main__':
+  main()

TestCases/parallel_regression.py

@@ -75,9 +75,17 @@ def main():
     cfd_flamelet_h2.cfg_file = "laminar_premixed_h2_flame_cfd.cfg"
     cfd_flamelet_h2.test_iter = 5
     cfd_flamelet_h2.test_vals = [-8.036958, -8.372668, -1.842800, -9.388446]
-    cfd_flamelet_h2.new_output = True
     test_list.append(cfd_flamelet_h2)
 
+    # Flame ignition methods
+    flame_init_methods = TestCase("flame_init_methods")
+    flame_init_methods.cfg_dir = "flamelet/08_flame_init_method"
+    flame_init_methods.cfg_file ="flamelet_initialization.py" 
+    flame_init_methods.command = TestCase.Command("mpirun -np 2", "python", "flamelet_initialization.py")
+    flame_init_methods.test_vals=[1600, 1600]
+    flame_init_methods.new_output = True 
+    test_list.append(flame_init_methods)
+
     #########################
     ## NEMO solver ###
     #########################

config_template.cfg

@@ -960,6 +960,16 @@ PREFERENTIAL_DIFFUSION= NO
 % (x8) = Thickness of the 'burnt' zone, after this length, the conditions will be 'unburnt' again.
 FLAME_INIT= (0.004, 0.0, 0.0, 1.0, 0.0, 0.0, 0.2e-3, 1.0)
 
+% Flame ignition temperature in Kelvin.
+% Specify the temperature on the 'burnt' side of the flame when using FLAME_FRONT initialization
+% or the maximum temperature within the spark region when using SPARK initialization.
+% The initial value of the progress variable on the 'burnt' side of the flame front is calculated through
+% inverse regression of the ignition temperature or when reaching the upper limit of the progress variable
+% in the flamelet data space.
+% The SPARK_REACTION_RATES are applied within the spark region when the local temperature is below the
+% flame ignition temperature.
+FLAME_INIT_IGNITION= 5000.0
+
 % SPARK initialization
 % 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

meson_scripts/init.py

@@ -74,7 +74,7 @@ def init_submodules(
     github_repo_mel = "https://github.com/pcarruscag/MEL"
     sha_version_fado = "ce7ee018e4e699af5028d69baa1939fea290e18a"
     github_repo_fado = "https://github.com/pcarruscag/FADO"
-    sha_version_mlpcpp = "02f2cb9dde791074858e11ac091f7c4df9c6af65"
+    sha_version_mlpcpp = "e23facf388902f262fbe7ba3bcc84d36c85350b9"
     github_repo_mlpcpp = "https://github.com/EvertBunschoten/MLPCpp"
     sha_version_eigen = "d71c30c47858effcbd39967097a2d99ee48db464"
     github_repo_eigen = "https://gitlab.com/libeigen/eigen.git"