[WIP] Fix flame initialization

Common/include/option_structure.hpp

@@ -1461,6 +1461,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;
 }

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