Flexible setup of flamelet boundary conditions

Common/include/CConfig.hpp

@@ -10265,4 +10265,12 @@ class CConfig {
    */
   const FluidFlamelet_ParsedOptions& GetFlameletParsedOptions() const { return flamelet_ParsedOptions; }
 
+  /*!
+   * \brief Get the enthalpy BC mode for the flamelet solver.
+   * FLOW_MARKERS: derive enthalpy BCs from MARKER_ISOTHERMAL/MARKER_HEATFLUX/MARKER_INLET (temperature-based).
+   * SPECIES_MARKERS: take enthalpy BCs directly from MARKER_WALL_SPECIES/MARKER_INLET_SPECIES.
+   * \return FLAMELET_ENTHALPY_BC enum value.
+   */
+  FLAMELET_ENTHALPY_BC GetFlamelet_Enthalpy_BC() const { return flamelet_ParsedOptions.enthalpy_bc; }
+
 };

Common/include/option_structure.hpp

@@ -1441,12 +1441,30 @@ static const MapType<std::string, FLAMELET_INIT_TYPE> Flamelet_Init_Map = {
   MakePair("SPARK", FLAMELET_INIT_TYPE::SPARK)
 };
 
+/*!
+ * \brief Selects the source of wall/inlet enthalpy boundary conditions for the flamelet solver.
+ * SPECIES_MARKERS (default): inlet H is taken directly from MARKER_INLET_SPECIES; wall enthalpy BC
+ * is obtained from MARKER_WALL_SPECIES.
+ * FLOW_MARKERS: inlet H is derived from the MARKER_INLET temperature via a Newton iteration on the
+ * LUT (reverse lookup using Z,T) from MARKER_ISOTHERMAL or MARKER_HEATFLUX.
+ */
+enum class FLAMELET_ENTHALPY_BC {
+  FLOW_MARKERS,    /*!< \brief Derive inlet H from MARKER_INLET T (LUT Newton); walls from MARKER_ISOTHERMAL/MARKER_HEATFLUX. */
+  SPECIES_MARKERS, /*!< \brief Take inlet H directly from MARKER_INLET_SPECIES; walls from MARKER_WALL_SPECIES (default). */
+};
+
+static const MapType<std::string, FLAMELET_ENTHALPY_BC> Flamelet_Enthalpy_BC_Map = {
+  MakePair("FLOW_MARKERS",    FLAMELET_ENTHALPY_BC::FLOW_MARKERS)
+  MakePair("SPECIES_MARKERS", FLAMELET_ENTHALPY_BC::SPECIES_MARKERS)
+};
+
 /*!
  * \brief Structure containing parsed options for flamelet fluid model.
  */
 struct FluidFlamelet_ParsedOptions {
   ///TODO: Add python wrapper initialization option
   FLAMELET_INIT_TYPE ignition_method = FLAMELET_INIT_TYPE::NONE; /*!< \brief Method for solution ignition for flamelet problems. */
+  FLAMELET_ENTHALPY_BC enthalpy_bc = FLAMELET_ENTHALPY_BC::SPECIES_MARKERS; /*!< \brief Source of enthalpy BCs: species markers (default, backward-compatible) or flow markers. */
   unsigned short n_scalars = 0;       /*!< \brief Number of transported scalars for flamelet LUT approach. */
   unsigned short n_lookups = 0;       /*!< \brief Number of lookup variables, for visualization only. */
   unsigned short n_table_sources = 0; /*!< \brief Number of transported scalar source terms for LUT. */

Common/src/CConfig.cpp

@@ -1410,6 +1410,8 @@ void CConfig::SetConfig_Options() {
 
   /*!\brief FLAME_INIT_METHOD \n DESCRIPTION: Ignition method for flamelet solver \n DEFAULT: no ignition; cold flow only. */
   addEnumOption("FLAME_INIT_METHOD", flamelet_ParsedOptions.ignition_method, Flamelet_Init_Map, FLAMELET_INIT_TYPE::NONE);
+  /*!\brief FLAME_ENTHALPY_BC \n DESCRIPTION: enthalpy BC from species markers (default, backward-compatible) or flow markers \n DEFAULT: SPECIES_MARKERS \ingroup Config */
+  addEnumOption("FLAME_ENTHALPY_BC", flamelet_ParsedOptions.enthalpy_bc, Flamelet_Enthalpy_BC_Map, FLAMELET_ENTHALPY_BC::SPECIES_MARKERS);
   /*!\brief FLAME_INIT \n DESCRIPTION: flame front initialization using the flamelet model \ingroup Config*/
   addDoubleArrayOption("FLAME_INIT", flamelet_ParsedOptions.flame_init.size(), false, flamelet_ParsedOptions.flame_init.begin());
 

SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp

@@ -162,6 +162,21 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
   void BC_Isothermal_Wall(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics,
                           CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) override;
 
+  /*!
+   * \brief Impose the heat-flux wall boundary condition on the flamelet enthalpy scalar.
+   * When FLAMELET_ENTHALPY_BC= FLOW_MARKERS (default), the heat flux value from MARKER_HEATFLUX
+   * is applied as a Neumann boundary condition on the enthalpy scalar H.
+   * When FLAMELET_ENTHALPY_BC= SPECIES_MARKERS, use flux/value from MARKER_WALL_SPECIES.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] conv_numerics - Description of the numerical method.
+   * \param[in] visc_numerics - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] val_marker - Surface marker where the boundary condition is applied.
+   */
+  void BC_HeatFlux_Wall(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics,
+                        CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) override;
+
   /*!
    * \brief Impose the inlet boundary condition.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp

@@ -38,6 +38,7 @@ class CSpeciesFlameletVariable final : public CSpeciesVariable {
   MatrixType source_scalar; /*!< \brief Vector of the source terms from the lookup table for each scalar equation */
   MatrixType lookup_scalar; /*!< \brief Vector of the source terms from the lookup table for each scalar equation */
   su2vector<unsigned short> table_misses; /*!< \brief Vector of lookup table misses. */
+  MatrixType source_cons_jac; /*!< \brief Consumption-rate Jacobian dS_aux_i/dY_aux_i = source_cons_i, one column per user scalar. */
 
  public:
   /*!
@@ -87,4 +88,18 @@ class CSpeciesFlameletVariable final : public CSpeciesVariable {
   inline void SetTableMisses(unsigned long iPoint, unsigned short misses) override { table_misses[iPoint] = misses; }
 
   inline unsigned short GetTableMisses(unsigned long iPoint) const override { return table_misses[iPoint]; }
+
+  /*!
+   * \brief Store the consumption-rate Jacobian dS_aux_i/dY_aux_i = source_cons_i for user scalar i_aux.
+   */
+  inline void SetAuxSourceCons(unsigned long iPoint, unsigned long i_aux, su2double val) {
+    source_cons_jac(iPoint, i_aux) = val;
+  }
+
+  /*!
+   * \brief Get the consumption-rate Jacobian entry for user scalar i_aux at iPoint.
+   */
+  inline su2double GetAuxSourceCons(unsigned long iPoint, unsigned long i_aux) const {
+    return source_cons_jac(iPoint, i_aux);
+  }
 };

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1183,7 +1183,7 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co
 
     if (LD2_Scheme) {
       numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(jPoint));
-      if (!geometry->nodes->GetPeriodicBoundary(iPoint) || (geometry->nodes->GetPeriodicBoundary(iPoint) 
+      if (!geometry->nodes->GetPeriodicBoundary(iPoint) || (geometry->nodes->GetPeriodicBoundary(iPoint)
           && !geometry->nodes->GetPeriodicBoundary(jPoint))) {
         numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint));
       } else {
@@ -2535,11 +2535,19 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
     if (species_model) scalar_inlet = config->GetInlet_SpeciesVal(config->GetMarker_All_TagBound(val_marker));
     CFluidModel* auxFluidModel = solver_container[FLOW_SOL]->GetFluidModel();
     auxFluidModel->SetTDState_T(V_inlet[prim_idx.Temperature()], scalar_inlet);
-    V_inlet[prim_idx.Enthalpy()] = auxFluidModel->GetEnthalpy();
+
+    /*--- For the flamelet model with FLOW_MARKERS enthalpy BC, we obtain the inlet enthalpy
+     from the flamelet species solver  With SPECIES_MARKERS, the enthalpy in MARKER_INLET_SPECIES
+     is used directly. ---*/
+    if (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET &&
+        config->GetFlamelet_Enthalpy_BC() == FLAMELET_ENTHALPY_BC::FLOW_MARKERS)
+      V_inlet[prim_idx.Enthalpy()] = nodes->GetEnthalpy(iPoint);
+    else
+      V_inlet[prim_idx.Enthalpy()] = auxFluidModel->GetEnthalpy();
 
     /*--- Access density at the node. This is either constant by
-      construction, or will be set fixed implicitly by the temperature
-      and equation of state. ---*/
+     construction, or will be set fixed implicitly by the temperature
+     and equation of state. ---*/
 
     V_inlet[prim_idx.Density()] = nodes->GetDensity(iPoint);
 

SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp

@@ -395,11 +395,28 @@ void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solv
                                              CNumerics** numerics_container, CConfig* config, unsigned short iMesh) {
   SU2_ZONE_SCOPED
 
+  const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  const auto n_CV = flamelet_config_options.n_control_vars;
+  const auto n_aux = flamelet_config_options.n_user_scalars;
+  const auto* fn = static_cast<const CSpeciesFlameletVariable*>(nodes);
+
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto i_point = 0u; i_point < nPointDomain; i_point++) {
+    const su2double volume = geometry->nodes->GetVolume(i_point);
+
     /*--- 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) -= nodes->GetScalarSources(i_point)[i_var] * volume;
+    }
+
+    /*--- Implicit: analytic Jacobian for auxiliary species from the split source form.
+     *   S_aux_i = source_prod_i + source_cons_i * Y_aux_i
+     *   dS_aux_i/dY_aux_i = source_cons_i
+     *   J_ii += -source_cons_i * V ---*/
+    if (implicit) {
+      for (auto i_aux = 0u; i_aux < n_aux; i_aux++) {
+        Jacobian.AddVal2Diag(i_point, n_CV + i_aux, -fn->GetAuxSourceCons(i_point, i_aux) * volume);
+      }
     }
   }
   END_SU2_OMP_FOR
@@ -409,17 +426,58 @@ void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solv
 
 }
 
+void CSpeciesFlameletSolver::BC_HeatFlux_Wall(CGeometry* geometry, CSolver** solver_container,
+                                               CNumerics* conv_numerics, CNumerics* visc_numerics,
+                                               CConfig* config, unsigned short val_marker) {
+  SU2_ZONE_SCOPED
+
+  /*--- In FLOW_MARKERS mode: read MARKER_HEATFLUX.
+   In SPECIES_MARKERS mode: read flux/value from MARKER_WALL_SPECIES. ---*/
+
+  if (config->GetFlamelet_Enthalpy_BC() != FLAMELET_ENTHALPY_BC::FLOW_MARKERS) {
+    CSpeciesSolver::BC_HeatFlux_Wall(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker);
+    return;
+  }
+
+  const string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
+  su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag);
+  if (config->GetIntegrated_HeatFlux())
+    Wall_HeatFlux /= geometry->GetSurfaceArea(config, val_marker);
+
+  SU2_OMP_FOR_DYN(OMP_MIN_SIZE)
+  for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) {
+    const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
+    if (!geometry->nodes->GetDomain(iPoint)) continue;
+
+    const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal();
+    const su2double Area = GeometryToolbox::Norm(nDim, Normal);
+
+    /*--- Neumann condition: q_wall is the prescribed heat flux (W/m^2, positive into domain).
+     This adds a source term dH/dn * lambda = q_wall to the enthalpy residual. ---*/
+    LinSysRes(iPoint, I_ENTH) -= Wall_HeatFlux * Area;
+  }
+  END_SU2_OMP_FOR
+}
+
 void CSpeciesFlameletSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics,
                                       CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) {
   SU2_ZONE_SCOPED
   string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
 
-  su2double temp_inlet = config->GetInletTtotal(Marker_Tag);
-
-  /*--- We compute inlet enthalpy from the temperature and progress variable. ---*/
   su2double enth_inlet;
-  GetEnthFromTemp(solver_container[FLOW_SOL]->GetFluidModel(), temp_inlet, config->GetInlet_SpeciesVal(Marker_Tag),
-                  &enth_inlet);
+  if (config->GetFlamelet_Enthalpy_BC() == FLAMELET_ENTHALPY_BC::FLOW_MARKERS) {
+    /*--- Derive inlet enthalpy from MARKER_INLET temperature via Newton iteration on the LUT.
+     This ensures the enthalpy is thermodynamically consistent with the prescribed temperature,
+     regardless of the value given in MARKER_INLET_SPECIES. ---*/
+    su2double temp_inlet = config->GetInletTtotal(Marker_Tag);
+    GetEnthFromTemp(solver_container[FLOW_SOL]->GetFluidModel(), temp_inlet,
+                    config->GetInlet_SpeciesVal(Marker_Tag), &enth_inlet);
+  } else {
+    /*--- Use the enthalpy value directly from MARKER_INLET_SPECIES (default).
+     The user is responsible for providing a thermodynamically consistent value. ---*/
+    enth_inlet = config->GetInlet_SpeciesVal(Marker_Tag)[I_ENTH];
+  }
+
   SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
   for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) {
     Inlet_SpeciesVars[val_marker][iVertex][I_ENTH] = enth_inlet;
@@ -552,6 +610,8 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF
     su2double source_prod = table_sources[flamelet_config_options.n_control_vars + 2 * i_aux];
     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;
+    /*--- Store the analytic Jacobian dS_aux/dY_aux = source_cons for implicit treatment. ---*/
+    static_cast<CSpeciesFlameletVariable*>(nodes)->SetAuxSourceCons(iPoint, i_aux, source_cons);
   }
   for (auto i_scalar = 0u; i_scalar < nVar; i_scalar++)
     nodes->SetScalarSource(iPoint, i_scalar, source_scalar[i_scalar]);

SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp

@@ -53,6 +53,7 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf,
   source_scalar.resize(nPoint, flamelet_config_options.n_scalars) = su2double(0.0);
   lookup_scalar.resize(nPoint, flamelet_config_options.n_lookups) = su2double(0.0);
   table_misses.resize(nPoint) = 0;
+  source_cons_jac.resize(nPoint, flamelet_config_options.n_user_scalars) = su2double(0.0);
 
   if (flamelet_config_options.preferential_diffusion) {
     AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS) = su2double(0.0);

config_template.cfg

@@ -1,6 +1,6 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %                                                                              %
-% SU2 configuration file                                                       %
+% SU2 configuration file config_template.cfg                                   %
 % Case description: _________________________________________________________  %
 % Author: ___________________________________________________________________  %
 % Institution: ______________________________________________________________  %
@@ -943,7 +943,12 @@ SPECIES_CLIPPING_MIN= 0.0
 % for definition of the flamelet manifold. Either LUT or MLP can be used as options.
 % If the terms "Beta_ProgVar", "Beta_Enth_Thermal", "Beta_Enth", and "Beta_Mixfrac" are
 % found in the variables list of the manifold, preferential diffusion is assumed.
+% The flamelet method uses the same boundary conditions as species transport.
 
+% FLOW_MARKERS: flamelet method uses temperature BC from MARKER_INLET and MARKER_HEATFLUX/MARKER_ISOTHERMAL.
+% It then does a reverse lookup for the enthalpy.
+% SPECIES_MARKERS: flamelet method uses enthalpy from MARKER_INLET_SPECIES and MARKER_WALL_SPECIES.
+FLAME_ENTHALPY_BC= SPECIES_MARKERS
 %
 % Names of the user defined (auxiliary) transport equations.
 USER_SCALAR_NAMES= (Y-CO)