Flexible setup of flamelet boundary conditions (#2814)

* flexible BC for flamelet equations
* add jacobian of the split source terms, rename FLAMELET->FLAME
* add python wrapper to flamelet enthalpy

---------

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>

Author: bigfooted

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,10 @@ 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 for thermal walls. FLOW_MARKERS (default): enthalpy derived
+   from MARKER_ISOTHERMAL temperature or MARKER_HEATFLUX via GetEnthFromTemp. SPECIES_MARKERS: enthalpy and all other
+   scalars taken directly from MARKER_WALL_SPECIES or the Python wrapper (SetMarkerCustomScalar). \n DEFAULT: FLOW_MARKERS \ingroup Config */
+  addEnumOption("FLAME_ENTHALPY_BC", flamelet_ParsedOptions.enthalpy_bc, Flamelet_Enthalpy_BC_Map, FLAMELET_ENTHALPY_BC::FLOW_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/flamelet_python_bc.md

@@ -0,0 +1,225 @@
+# Flamelet Solver — Python Wrapper Wall Boundary Conditions
+
+This document explains how to use the Python wrapper (`pysu2`) to set
+custom per-vertex boundary conditions for the flamelet (FGM) scalar transport
+solver, and how the setup differs between the two `FLAMELET_ENTHALPY_BC` modes.
+
+---
+
+## Background
+
+The flamelet solver transports `nVar` scalar variables. Their order is:
+
+| Index | Variable |
+|-------|----------|
+| 0 | `PROGVAR` (progress variable) |
+| 1 | `ENTH` (total enthalpy) |
+| 2 … n\_CV−1 | additional control variables (e.g. `MIXFRAC`) |
+| n\_CV … nVar−1 | auxiliary / user-defined species |
+
+Two modes control how the enthalpy boundary condition is derived at thermal
+walls (`MARKER_ISOTHERMAL`, `MARKER_HEATFLUX`):
+
+- **`FLOW_MARKERS`** — enthalpy is derived from the flow thermal field (wall
+  temperature via `MARKER_ISOTHERMAL`, or wall heat flux via `MARKER_HEATFLUX`).
+- **`SPECIES_MARKERS`** (default) — enthalpy and all other scalars are taken
+  directly from `MARKER_WALL_SPECIES`, or from the Python wrapper when
+  `MARKER_PYTHON_CUSTOM` is active.
+
+---
+
+## Mode 1: `FLAMELET_ENTHALPY_BC = FLOW_MARKERS`
+
+### How it works
+
+The C++ function `BC_HeatFlux_Wall` reads the wall heat flux and applies it as
+a Neumann condition on `I_ENTH`.  When the marker is also listed in
+`MARKER_PYTHON_CUSTOM`, it instead reads a **per-vertex heat flux** set by the
+Python wrapper.
+
+```
+driver.SetMarkerCustomNormalHeatFlux(iMarker, iVertex, q_wall)
+    → geometry->CustomBoundaryHeatFlux[iMarker][iVertex]
+        → BC_HeatFlux_Wall reads geometry->GetCustomBoundaryHeatFlux(...)
+```
+
+> `BC_Isothermal_Wall` in `FLOW_MARKERS` mode always converts the wall
+> temperature (from `MARKER_ISOTHERMAL`) to enthalpy via `GetEnthFromTemp`.
+> Per-vertex customisation of the enthalpy Dirichlet value is not supported in
+> this mode — use `SPECIES_MARKERS` instead.
+
+### Limitations in this mode
+
+- Only `I_ENTH` is reachable via the Python wrapper.
+- Non-enthalpy scalars (`PROGVAR`, auxiliary species) receive implicit zero
+  Neumann at the wall and cannot be set via Python.
+
+### config.cfg
+
+```cfg
+FLAMELET_ENTHALPY_BC = FLOW_MARKERS
+
+% Flow solver wall BC — sets KindBC and provides the fallback heat flux value
+% used when MARKER_PYTHON_CUSTOM is not active.
+MARKER_HEATFLUX = ( burner_wall, 0.0 )
+
+% Enable the per-vertex Python override path.
+MARKER_PYTHON_CUSTOM = ( burner_wall )
+```
+
+### run.py
+
+```python
+import pysu2
+
+driver = pysu2.CSinglezoneDriver("config.cfg", 1, False)
+
+marker_name = "burner_wall"
+iMarker     = list(driver.GetMarkerTags()).index(marker_name)
+n_vertex    = driver.GetNumberMarkerNodes(iMarker)
+
+def compute_heat_flux(coord):
+    """Heat flux [W/m²], positive = into the domain."""
+    import math
+    r = math.sqrt(coord[0]**2 + coord[1]**2)
+    return -5000.0 * math.exp(-r**2 / 0.01)   # Gaussian profile, heat out
+
+for iteration in range(driver.GetNumberIter()):
+
+    for iVertex in range(n_vertex):
+        node_id = driver.GetMarkerNode(iMarker, iVertex)
+        coord   = driver.GetInitialMeshCoord(node_id)   # [x, y(, z)]
+        q_wall  = compute_heat_flux(coord)
+        driver.SetMarkerCustomNormalHeatFlux(iMarker, iVertex, q_wall)
+
+    driver.Preprocess(iteration)
+    driver.Run()
+    driver.Postprocess()
+    driver.Monitor(iteration)
+    driver.Output(iteration)
+
+driver.Finalize()
+```
+
+---
+
+## Mode 2: `FLAMELET_ENTHALPY_BC = SPECIES_MARKERS` (default)
+
+### How it works
+
+Both `BC_HeatFlux_Wall` and `BC_Isothermal_Wall` delegate immediately to
+`CSpeciesSolver::BC_Wall_Generic`, which processes **all `nVar` scalars**
+independently.  When `MARKER_PYTHON_CUSTOM` is active, the value for each
+scalar is taken from the array set by `driver.SetMarkerCustomScalar()`.
+
+```
+driver.SetMarkerCustomScalar(iMarker, iVertex, [val_0, val_1, ..., val_nVar-1])
+    → CustomBoundaryScalar[iMarker](iVertex, iVar)
+        → BC_Wall_Generic reads CustomBoundaryScalar for every iVar
+```
+
+The **type** of BC for each scalar (Dirichlet `VALUE` or Neumann `FLUX`) is
+read from `MARKER_WALL_SPECIES` in the config and cannot be changed from
+Python.  The Python wrapper only overrides the **magnitude** per vertex.
+
+### config.cfg
+
+The example below uses `nVar = 3` (PROGVAR, ENTH, one auxiliary species).
+Adjust the number of `BC_TYPE, value` pairs to match the actual number of
+scalar variables in your setup.
+
+```cfg
+FLAMELET_ENTHALPY_BC = SPECIES_MARKERS
+
+% Flow solver wall BC — still required to set KindBC.
+% The choice between MARKER_ISOTHERMAL and MARKER_HEATFLUX only affects
+% the flow solver; the species solver uses BC_Wall_Generic in both cases.
+MARKER_ISOTHERMAL = ( burner_wall, 300.0 )
+
+% Per-variable BC types for the species solver.
+% Format: (marker_name, BC_TYPE, fallback_value, BC_TYPE, fallback_value, ...)
+%         One BC_TYPE+value pair per scalar variable, in index order.
+% BC_TYPE: FLUX  → Neumann (value is normal flux density [unit/m²])
+%          VALUE → Dirichlet (value is the scalar value at the wall)
+%
+% When MARKER_PYTHON_CUSTOM is active, the fallback_value is ignored and
+% the Python-supplied value is used instead.  The BC_TYPE is always from config.
+%
+% idx 0 PROGVAR : zero Neumann  (no progress variable source at the wall)
+% idx 1 ENTH    : Dirichlet     (enthalpy value set per-vertex by Python)
+% idx 2 aux_1   : zero Neumann  (no auxiliary species flux at the wall)
+MARKER_WALL_SPECIES = ( burner_wall, FLUX, 0.0, VALUE, 0.0, FLUX, 0.0 )
+
+% Enable the per-vertex Python override path.
+MARKER_PYTHON_CUSTOM = ( burner_wall )
+```
+
+### run.py
+
+```python
+import pysu2
+
+driver = pysu2.CSinglezoneDriver("config.cfg", 1, False)
+
+marker_name = "burner_wall"
+iMarker     = list(driver.GetMarkerTags()).index(marker_name)
+n_vertex    = driver.GetNumberMarkerNodes(iMarker)
+
+def compute_wall_enthalpy(coord):
+    """Return wall enthalpy [J/kg] at this vertex."""
+    # Example: uniform value corresponding to T=300 K from your LUT.
+    return 3.5e5
+
+for iteration in range(driver.GetNumberIter()):
+
+    for iVertex in range(n_vertex):
+        node_id = driver.GetMarkerNode(iMarker, iVertex)
+        coord   = driver.GetInitialMeshCoord(node_id)
+
+        enth_wall = compute_wall_enthalpy(coord)
+
+        # Pass one value per scalar variable (length = nVar).
+        # PROGVAR (idx 0): 0.0 → applied as FLUX → zero Neumann (BC_TYPE from config)
+        # ENTH    (idx 1): enth_wall → applied as VALUE → Dirichlet
+        # aux_1   (idx 2): 0.0 → applied as FLUX → zero Neumann
+        driver.SetMarkerCustomScalar(iMarker, iVertex, [0.0, enth_wall, 0.0])
+
+    driver.Preprocess(iteration)
+    driver.Run()
+    driver.Postprocess()
+    driver.Monitor(iteration)
+    driver.Output(iteration)
+
+driver.Finalize()
+```
+
+---
+
+## Comparison
+
+| | `FLOW_MARKERS` | `SPECIES_MARKERS` |
+|---|---|---|
+| **Python API** | `SetMarkerCustomNormalHeatFlux` | `SetMarkerCustomScalar` |
+| **Storage** | `geometry->CustomBoundaryHeatFlux` | `CustomBoundaryScalar[iMarker](iVertex, iVar)` |
+| **What you set** | Per-vertex heat flux (W/m²) for `I_ENTH` only | Per-vertex values for **all** `nVar` scalars |
+| **BC type control** | Hard-coded Neumann in C++ | Per-variable via `MARKER_WALL_SPECIES` in config |
+| **`MARKER_ISOTHERMAL` walls** | Enthalpy computed from T via `GetEnthFromTemp`; not overridable by Python | Enthalpy set as Dirichlet `VALUE` from Python |
+| **Non-enthalpy scalars** | Not reachable via Python | All controlled via the same `SetMarkerCustomScalar` call |
+| **Required config keys** | `MARKER_HEATFLUX`, `MARKER_PYTHON_CUSTOM` | `MARKER_ISOTHERMAL` or `MARKER_HEATFLUX`, `MARKER_WALL_SPECIES`, `MARKER_PYTHON_CUSTOM` |
+
+---
+
+## Notes
+
+- `MARKER_PYTHON_CUSTOM` is a purely additive flag. A marker can appear in
+  both `MARKER_HEATFLUX` (or `MARKER_ISOTHERMAL`) and `MARKER_PYTHON_CUSTOM`
+  simultaneously.
+- The fallback values in `MARKER_WALL_SPECIES` are used when `py_custom` is
+  false (e.g., during a plain SU2\_CFD run without the Python driver). They
+  allow the same config to be used both ways.
+- In `SPECIES_MARKERS` mode, `SetMarkerCustomScalar` must supply a vector of
+  exactly `nVar` values. Indices that correspond to `FLUX`-type variables
+  should receive `0.0` unless you intentionally want a non-zero flux there.
+- The `CHT` (conjugate heat transfer) path (`BC_ConjugateHeat_Interface`)
+  always uses `BC_Isothermal_Wall_Generic` and is unaffected by
+  `FLAMELET_ENTHALPY_BC`. Python custom BCs are not applicable to CHT markers.

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);