Source term method for preferential diffusion

Common/include/containers/CLookUpTable.hpp

@@ -82,6 +82,10 @@ class CLookUpTable {
 
   double memory_footprint_data = 0; /*!< \brief Memory footprint of the loaded table data. */
 
+  su2double hull_miss_dist_ = 0.0;  /*!< \brief Worst-case normalised distance to hull boundary (used to select the reporting level). */
+  su2double cv1_hull_dev_  = 0.0;  /*!< \brief Signed physical deviation in CV1 (query minus nearest hull node) at the worst-miss level. */
+  su2double cv2_hull_dev_  = 0.0;  /*!< \brief Signed physical deviation in CV2 (query minus nearest hull node) at the worst-miss level. */
+
   /*! \brief
    * Holds all connectivity data stored in the table for each level. First index
    * addresses the variable while second index addresses the point.
@@ -370,6 +374,18 @@ class CLookUpTable {
    */
   std::pair<unsigned long, unsigned long> FindInclusionLevels(const su2double val_CV3);
 
+  /*!
+   * \brief Distance from val_CV3 to the nearest table Z level (in physical Z units).
+   *        Returns 0 for 2D tables or when only one Z level exists.
+   */
+  inline su2double GetDistanceToNearestZLevel(su2double val_CV3) const {
+    if (table_dim < 3 || n_table_levels < 2) return 0.0;
+    auto it = std::lower_bound(z_values_levels.begin(), z_values_levels.end(), val_CV3);
+    su2double d_hi = (it != z_values_levels.end())   ? abs(*it - val_CV3) : su2double(1e30);
+    su2double d_lo = (it != z_values_levels.begin()) ? abs(*std::prev(it) - val_CV3) : su2double(1e30);
+    return (d_lo < d_hi) ? d_lo : d_hi;
+  }
+
   /*!
    * \brief Determine the minimum and maximum value of the second controlling variable.
    * \returns Pair of minimum and maximum value of controlling variable 2.
@@ -386,6 +402,37 @@ class CLookUpTable {
     return limits_table_x[i_level];
   }
 
+  /*!
+   * \brief Get the global Z-dimension limits of the table (min, max) as values.
+   * \returns Pair of (z_min, z_max); both zero for 2D tables.
+   */
+  inline std::pair<su2double, su2double> GetTableLimitsZ() const {
+    if (table_dim < 3) return {su2double(0), su2double(0)};
+    return {*limits_table_z.first, *limits_table_z.second};
+  }
+
+  /*!
+   * \brief Get the number of table levels (Z dimension slices).
+   */
+  inline unsigned long GetNTableLevels() const { return n_table_levels; }
+
+  /*!
+   * \brief Reset the hull-miss accumulators (call once per point before all lookups).
+   */
+  void ResetHullMissDistance() { hull_miss_dist_ = 0.0; cv1_hull_dev_ = 0.0; cv2_hull_dev_ = 0.0; }
+
+  /*!
+   * \brief Signed physical deviation in CV1 (query minus nearest hull node) at the worst-miss Z level.
+   *        Returns 0 if the last query was inside the hull.
+   */
+  su2double GetHullMissCV1Dev() const { return cv1_hull_dev_; }
+
+  /*!
+   * \brief Signed physical deviation in CV2 (query minus nearest hull node) at the worst-miss Z level.
+   *        Returns 0 if the last query was inside the hull.
+   */
+  su2double GetHullMissCV2Dev() const { return cv2_hull_dev_; }
+
   /*!
    * \brief Check whether requested set of variables are included in the table.
    */

Common/include/option_structure.hpp

@@ -1546,6 +1546,23 @@ enum FLAMELET_PREF_DIFF_SCALARS {
   N_BETA_TERMS,         /*!< \brief Total number of preferential diffusion scalars. */
 };
 
+/*!
+ * \brief Preferential diffusion flux term index within one control variable's block of
+ * coefficients: the molecular coefficients of the major species occupy [0, n_major_species),
+ * and the thermal (Soret) coefficient occupies index n_major_species. The major species
+ * themselves are configured at run time with PREFERENTIAL_DIFFUSION_MAJOR_SPECIES, so their
+ * count is not a compile-time constant.
+ */
+static inline unsigned short FlameletPDThermalTerm(unsigned short n_major_species) { return n_major_species; }
+
+/*!
+ * \brief Number of preferential diffusion flux coefficients per control variable: one molecular
+ * coefficient per major species plus the single thermal (Soret) coefficient.
+ */
+static inline unsigned short FlameletPDTermsPerCV(unsigned short n_major_species) {
+  return n_major_species + 1;
+}
+
 /*!
  * \brief Flame initialization options for the flamelet solver.
  */
@@ -1578,6 +1595,24 @@ static const MapType<std::string, FLAMELET_ENTHALPY_BC> Flamelet_Enthalpy_BC_Map
   MakePair("SPECIES_MARKERS", FLAMELET_ENTHALPY_BC::SPECIES_MARKERS)
 };
 
+/*!
+ * \brief Selects which preferential diffusion method is active in the flamelet scalar solver.
+ * BETA_CORRECTION (default): β-scalar viscous flux correction (constant-Lewis formulation,
+ * Mukundakumar et al.).
+ * SOURCE_TERM: major-species model B2 of Schepers & van Oijen, C&F 280 (2025) 114332 —
+ * runtime Eq. (14) fluxes (molecular D_{phi_k,i} grad(Y_i) for i in {H2, H2O, H} plus thermal
+ * D^T_{phi_k} grad(T)) combined with the Eq. (16) closure source of the non-major species.
+ */
+enum class FLAMELET_PD_METHOD {
+  BETA_CORRECTION, /*!< \brief β-scalar viscous flux correction applied to the diffusion operator. */
+  SOURCE_TERM,     /*!< \brief Major-species PD fluxes + non-major closure source terms. */
+};
+
+static const MapType<std::string, FLAMELET_PD_METHOD> Flamelet_PD_Method_Map = {
+  MakePair("BETA_CORRECTION", FLAMELET_PD_METHOD::BETA_CORRECTION)
+  MakePair("SOURCE_TERM",     FLAMELET_PD_METHOD::SOURCE_TERM)
+};
+
 /*!
  * \brief Structure containing parsed options for flamelet fluid model.
  */
@@ -1606,6 +1641,11 @@ struct FluidFlamelet_ParsedOptions {
   unsigned short nspark;           /*!< \brief Number of source terms for spark initialization. */
   bool preferential_diffusion = false;  /*!< \brief Preferential diffusion physics for flamelet solver.*/
   bool thickenedflame_correction{true}; /*!< \brief Thickened flame correction. */
+  unsigned short n_pd_major_species = 0;   /*!< \brief Number of preferential diffusion major species (SOURCE_TERM method). */
+  std::string* pd_major_species_names;     /*!< \brief Names of the preferential diffusion major species; the manifold
+                                                variables "D_<cv>_<species>" are composed from these. */
+  bool verbose_misses = false;                   /*!< \brief Print per-CV breakdown of manifold miss counts. */
+  FLAMELET_PD_METHOD pd_method = FLAMELET_PD_METHOD::BETA_CORRECTION; /*!< \brief Active PD term selection. */
   su2double Flame_T_ignition = 5000;    /*!< \brief Ignition temperature for the flame, used for initialization. */
 
 };

Common/src/CConfig.cpp

@@ -2394,6 +2394,20 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Enable preferential diffusion for FGM simulations. \n DEFAULT: false */
   addBoolOption("PREFERENTIAL_DIFFUSION", flamelet_ParsedOptions.preferential_diffusion, false);
 
+  /* DESCRIPTION: Print per-CV breakdown of manifold miss counts (Progress variable / Enthalpy / Mixture fraction / Hull). \n DEFAULT: false */
+  addBoolOption("FLAMELET_VERBOSE_MISSES", flamelet_ParsedOptions.verbose_misses, false);
+
+  /* DESCRIPTION: Active preferential diffusion terms. BETA_CORRECTION (default): β-scalar viscous flux
+   * correction. SOURCE_TERM: PD closure source terms for non-major species. COMBINED: both. */
+  addEnumOption("PREFERENTIAL_DIFFUSION_METHOD", flamelet_ParsedOptions.pd_method, Flamelet_PD_Method_Map, FLAMELET_PD_METHOD::BETA_CORRECTION);
+
+  /* DESCRIPTION: Major species carrying the resolved Eq. (14) preferential diffusion fluxes of the
+   * SOURCE_TERM method. The manifold variables are composed from these names, so the order and
+   * spelling must match the table: species "H2" implies "D_<cv>_H2" for every controlling variable
+   * <cv> and the mass fraction "Y-H2". Required when PREFERENTIAL_DIFFUSION_METHOD= SOURCE_TERM. */
+  addStringListOption("PREFERENTIAL_DIFFUSION_MAJOR_SPECIES", flamelet_ParsedOptions.n_pd_major_species,
+                      flamelet_ParsedOptions.pd_major_species_names);
+
   /*!\brief CONV_FILENAME \n DESCRIPTION: Output file convergence history (w/o extension) \n DEFAULT: history \ingroup Config*/
   addStringOption("CONV_FILENAME", Conv_FileName, string("history"));
   /*!\brief BREAKDOWN_FILENAME \n DESCRIPTION: Output file forces breakdown \ingroup Config*/
@@ -6019,6 +6033,18 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     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);
     }
+
+    /*--- The SOURCE_TERM preferential diffusion method resolves the Eq. (14) flux of a set of major
+     species at run time, so that set has to be named: the manifold variables looked up
+     ("D_<cv>_<species>", "Y-<species>") are composed from these names. Without them there is nothing
+     to look up and the method degenerates silently to no preferential diffusion at all. ---*/
+    if (flamelet_ParsedOptions.preferential_diffusion &&
+        flamelet_ParsedOptions.pd_method == FLAMELET_PD_METHOD::SOURCE_TERM &&
+        flamelet_ParsedOptions.n_pd_major_species == 0) {
+      SU2_MPI::Error("PREFERENTIAL_DIFFUSION_METHOD= SOURCE_TERM requires PREFERENTIAL_DIFFUSION_MAJOR_SPECIES "
+                     "to list the major species carried by the resolved preferential diffusion flux "
+                     "(e.g. PREFERENTIAL_DIFFUSION_MAJOR_SPECIES= (H2, H2O, H)).", CURRENT_FUNCTION);
+    }
   }
 
 }

SU2_CFD/include/fluid/CFluidFlamelet.hpp

@@ -54,6 +54,9 @@ class CFluidFlamelet final : public CFluidModel {
   unsigned short n_scalars, n_lookups, n_user_scalars, /*!< \brief number of passive reactant species. */
       n_control_vars;                                  /*!< \brief number of controlling variables. */
 
+  unsigned short n_pd_major_species = 0; /*!< \brief number of preferential diffusion major species carrying the
+                                              resolved Eq. (14) flux (SOURCE_TERM method, zero otherwise). */
+
   unsigned long extrapolation;
 
   INC_DENSITYMODEL density_model;
@@ -167,4 +170,23 @@ class CFluidFlamelet final : public CFluidModel {
    * \return Inclusion of preferential diffusion model.
    */
   inline bool GetPreferentialDiffusion() const override { return preferential_diffusion; }
+
+  /*!
+   * \brief Get the global bounds of all controlling variables over all table levels.
+   * Used for per-CV miss classification when FLAMELET_VERBOSE_MISSES is enabled.
+   */
+  void GetTableCVBounds(su2double& cv1_min, su2double& cv1_max,
+                        su2double& cv2_min, su2double& cv2_max,
+                        su2double& cv3_min, su2double& cv3_max) const;
+
+  void ResetHullMissDistance() override;
+  su2double GetHullMissCV1Dev() const override;
+  su2double GetHullMissCV2Dev() const override;
+
+  /*!
+   * \brief Distance from val_CV3 to the nearest table Z level (physical Z units).
+   *        Returns 0 for 2D tables.
+   */
+  su2double GetDistanceToNearestZLevel(su2double val_CV3) const;
+
 };

SU2_CFD/include/fluid/CFluidModel.hpp

@@ -398,6 +398,24 @@ class CFluidModel {
    */
   virtual bool GetPreferentialDiffusion() const { return false; }
 
+  /*!
+   * \brief Reset the hull-miss accumulators before the lookups for a new point.
+   *        No-op for fluid models without a convex-hull-based table.
+   */
+  virtual void ResetHullMissDistance() {}
+
+  /*!
+   * \brief Signed physical deviation in CV1 (query minus nearest hull node) at the worst-miss Z level.
+   *        Returns 0 for in-hull queries or fluid models without hull-based tables.
+   */
+  virtual su2double GetHullMissCV1Dev() const { return 0.0; }
+
+  /*!
+   * \brief Signed physical deviation in CV2 (query minus nearest hull node) at the worst-miss Z level.
+   *        Returns 0 for in-hull queries or fluid models without hull-based tables.
+   */
+  virtual su2double GetHullMissCV2Dev() const { return 0.0; }
+
   /*!
    * \brief Get number of Newton solver iterations.
    * \return Newton solver iteration count at termination.