diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index e6bf86d074d2..8299093e588d 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/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 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 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. */ diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 351e534671ae..12bb16794119 100644 --- a/Common/include/option_structure.hpp +++ b/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 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 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__" 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. */ }; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 1ab7ced2d354..9a344199d97f 100644 --- a/Common/src/CConfig.cpp +++ b/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__H2" for every controlling variable + * 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__", "Y-") 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); + } } } diff --git a/SU2_CFD/include/fluid/CFluidFlamelet.hpp b/SU2_CFD/include/fluid/CFluidFlamelet.hpp index 7b3cbeefa58e..6bec8f02f523 100644 --- a/SU2_CFD/include/fluid/CFluidFlamelet.hpp +++ b/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; + }; diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index 0e0c2e624bd0..e83b9b69835c 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/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. diff --git a/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp b/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp index dc634a235ca1..9861550451c0 100644 --- a/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp +++ b/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp @@ -28,6 +28,7 @@ #pragma once #include "species_edge_flux.hpp" +#include "../../variables/CSpeciesFlameletVariable.hpp" /*! * \class CScalarFlux_Flamelet @@ -50,7 +51,10 @@ class CScalarFlux_Flamelet final explicit CScalarFlux_Flamelet(const CConfig& config) : Base(config), preferentialDiffusion(config.GetFlameletParsedOptions().preferential_diffusion), - nControlVars(config.GetFlameletParsedOptions().n_control_vars) {} + nControlVars(config.GetFlameletParsedOptions().n_control_vars), + pdMethod(config.GetFlameletParsedOptions().pd_method), + nMajorSpecies(config.GetFlameletParsedOptions().n_pd_major_species), + pdTermsPerCV(FlameletPDTermsPerCV(config.GetFlameletParsedOptions().n_pd_major_species)) {} /*! * \brief Preferential diffusion, two terms with the shape of the ordinary diffusion but of @@ -68,6 +72,11 @@ class CScalarFlux_Flamelet final EdgeResidual& res) const { if (!preferentialDiffusion) return; + if (pdMethod == FLAMELET_PD_METHOD::SOURCE_TERM) { + sourceTermFluxes(idx, opt, iPoint, side_i, jPoint, side_j, rho, normal, vector_ij, res); + return; + } + const Double dist2_ij = fmax(squaredNorm(vector_ij), EPS); const Double proj_vector_ij = dot(vector_ij, normal) / dist2_ij; const Double proj_on_rho_i = proj_vector_ij / rho.i; @@ -137,6 +146,126 @@ class CScalarFlux_Flamelet final private: const bool preferentialDiffusion; const unsigned short nControlVars; + const FLAMELET_PD_METHOD pdMethod; + const unsigned short nMajorSpecies; + const unsigned short pdTermsPerCV; + + /*! + * \brief Resolved preferential diffusion fluxes of the SOURCE_TERM method, Eq. (14) of + * Schepers & van Oijen (C&F 280, 2025): for each controlling variable phi_k, + * J_{phi_k} = - sum_sp D_{phi_k,sp} grad(Y_sp) - (D^T_{phi_k}/T) grad(T), + * entering the transport equation as -div(J). The major species mass fractions are the + * solver's auxiliary variables, so their CFD-resolved gradients drive the flux, which is + * the point of tabulating a coefficient per species instead of pre-contracting them + * against the one-dimensional flamelet gradients. + * \note The tabulated coefficients are already density weighted (kg/(m s), or W/m for the + * enthalpy row), unlike the kinematic diffusivity the ordinary term uses, so they are + * averaged directly and NOT multiplied by rho. The 1/T of the thermal flux is not + * tabulated and is applied here with the local temperature. + * \note Explicit: the coefficients and the major species are manifold lookups, so there is no + * exact Jacobian with respect to the transported scalars. The implicit part is a + * stabilising self-diffusion added to the Jacobian alone, below. + */ + template + FORCEINLINE void sourceTermFluxes(const FlowIndices& idx, const ScalarFluxOptions& opt, Int iPoint, + const EdgeSide& side_i, Int jPoint, + const EdgeSide& side_j, const CPair& rho, + const Vector& normal, const Vector& vector_ij, + EdgeResidual& res) const { + const Double dist2_ij = fmax(squaredNorm(vector_ij), EPS); + const Double proj_vector_ij = dot(vector_ij, normal) / dist2_ij; + + const Double T_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.Temperature()); + const Double T_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.Temperature()); + const auto gradT_i = gatherVariables(iPoint, side_i.flowNodes->GetGradient_Primitive(), idx.Temperature()); + const auto gradT_j = gatherVariables(jPoint, side_j.flowNodes->GetGradient_Primitive(), idx.Temperature()); + + /*--- The scalar solvers are templated on CSpeciesVariable, which does not carry the Eq. (14) + coefficients. Only the flamelet solver instantiates this flux, and only on its own interior + edges, so its container is what these sides hold. ---*/ + const auto& pdCoeff_i = static_cast(side_i.scalarNodes).GetPDFluxCoeffs(); + const auto& pdCoeff_j = static_cast(side_j.scalarNodes).GetPDFluxCoeffs(); + + for (auto iScalar = 0u; iScalar < nControlVars; ++iScalar) { + /*--- Magnitude of the Eq. (14) flux on each side, accumulated as the terms are applied, for + the stabilising diffusivity below. ---*/ + Double fluxMag_i = 0.0, fluxMag_j = 0.0; + + /*--- Molecular terms, one per major species. ---*/ + for (auto iSp = 0u; iSp < nMajorSpecies; ++iSp) { + const auto col = iScalar * pdTermsPerCV + iSp; + + const Double Y_i = gatherVariables(iPoint, side_i.scalarNodes.GetAuxVar(), iSp); + const Double Y_j = gatherVariables(jPoint, side_j.scalarNodes.GetAuxVar(), iSp); + const auto gradY_i = gatherVariables(iPoint, side_i.scalarNodes.GetAuxVarGradient(), iSp); + const auto gradY_j = gatherVariables(jPoint, side_j.scalarNodes.GetAuxVarGradient(), iSp); + + const Double c_i = gatherVariables(iPoint, pdCoeff_i, col); + const Double c_j = gatherVariables(jPoint, pdCoeff_j, col); + const Double D = 0.5 * (c_i + c_j); + + const Double projGrad = projectedGradient(opt, gradY_i, gradY_j, Y_i, Y_j, normal, vector_ij, dist2_ij); + + res.flux_i(iScalar) -= D * projGrad; + if (!opt.oneSided) res.flux_j(iScalar) += D * projGrad; + + if (opt.implicit) { + fluxMag_i += fabs(c_i) * sqrt(fmax(squaredNorm(gradY_i), EPS)); + fluxMag_j += fabs(c_j) * sqrt(fmax(squaredNorm(gradY_j), EPS)); + } + } + + /*--- Thermal (Soret) term. ---*/ + const auto colT = iScalar * pdTermsPerCV + FlameletPDThermalTerm(nMajorSpecies); + const Double Dth = 0.5 * (gatherVariables(iPoint, pdCoeff_i, colT) / T_i + + gatherVariables(jPoint, pdCoeff_j, colT) / T_j); + + const Double projGradT = projectedGradient(opt, gradT_i, gradT_j, T_i, T_j, normal, vector_ij, dist2_ij); + + res.flux_i(iScalar) -= Dth * projGradT; + if (!opt.oneSided) res.flux_j(iScalar) += Dth * projGradT; + + /*--- Stabilisation. The fluxes above are explicit, and explicit diffusion carries a limit on + the pseudo time step that the flame front cells reach first. A self-diffusion of the + controlling variable is added to the Jacobian ALONE: the residual still holds exactly the + Eq. (14) fluxes, so the converged solution is untouched, while the positive definite + operator lifts that limit the way an implicit treatment of a real self-diffusion would. + D_stab is the baseline diffusivity scaled by how much larger the Eq. (14) coefficients are, + capped so a vanishing baseline cannot make it unbounded. ---*/ + if (opt.implicit) { + const Double Dbase_i = rho.i * gatherVariables(iPoint, side_i.scalarNodes.GetDiffusivity(), iScalar); + const Double Dbase_j = rho.j * gatherVariables(jPoint, side_j.scalarNodes.GetDiffusivity(), iScalar); + + /*--- Recast the flux bound as a self diffusivity, |J| <= D_stab |grad phi|, so the implicit + operator dominates the explicit flux without damping more than the flux itself warrants. + Scaling by the control variable's own gradient is what keeps it that tight: dropping it + leaves D_stab at the raw coefficient magnitude, which over damps and cancels the very CFL + gain the term exists to provide. Where grad(phi) vanishes while the cross gradients do not + the bound degenerates, so it is capped against the baseline diffusion. ---*/ + fluxMag_i += fabs(gatherVariables(iPoint, pdCoeff_i, colT) / T_i) * sqrt(fmax(squaredNorm(gradT_i), EPS)); + fluxMag_j += fabs(gatherVariables(jPoint, pdCoeff_j, colT) / T_j) * sqrt(fmax(squaredNorm(gradT_j), EPS)); + + const auto gradPhi_i = gatherVariables(iPoint, side_i.scalarNodes.GetGradient(), iScalar); + const auto gradPhi_j = gatherVariables(jPoint, side_j.scalarNodes.GetGradient(), iScalar); + + const Double Dstab_i = fmin(fluxMag_i / sqrt(fmax(squaredNorm(gradPhi_i), EPS)), C_STAB_MAX * Dbase_i); + const Double Dstab_j = fmin(fluxMag_j / sqrt(fmax(squaredNorm(gradPhi_j), EPS)), C_STAB_MAX * Dbase_j); + const Double Dstab = 0.5 * (Dstab_i + Dstab_j); + + res.jac_ii(iScalar, iScalar) += Dstab * proj_vector_ij / rho.i; + if (!opt.oneSided) { + res.jac_ij(iScalar, iScalar) -= Dstab * proj_vector_ij / rho.j; + res.jac_ji(iScalar, iScalar) -= Dstab * proj_vector_ij / rho.i; + res.jac_jj(iScalar, iScalar) += Dstab * proj_vector_ij / rho.j; + } + } + } + } + + /*! + * \brief Cap on the stabilising diffusivity, as a multiple of the baseline diffusion. + */ + static constexpr passivedouble C_STAB_MAX = 100.0; /*! * \brief Auxiliary variable holding the beta scalar of a controlling variable. diff --git a/SU2_CFD/include/output/CFlowOutput.hpp b/SU2_CFD/include/output/CFlowOutput.hpp index 7e6264f18476..da7c4d476b6e 100644 --- a/SU2_CFD/include/output/CFlowOutput.hpp +++ b/SU2_CFD/include/output/CFlowOutput.hpp @@ -37,6 +37,7 @@ struct CPrimitiveIndices; class CFlowOutput : public CFVMOutput{ protected: unsigned long lastInnerIter; + su2double flamelet_pv_range = 1.0; /*!< \brief Global PV range for C+ flame resolution index; refreshed each output write. */ /*! * \brief Constructor of the class @@ -168,6 +169,21 @@ class CFlowOutput : public CFVMOutput{ */ void SetVolumeOutputFieldsScalarMisc(const CConfig* config); + /*! + * \brief Compute the per-write global quantities the volume output needs (flamelet progress + * variable range for the C+ index). Collective; must not be reached per point. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - The container holding all solution data. + */ + void PrepareVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver) override; + + /*! + * \brief Add flamelet mesh quality diagnostic fields (FLAME_QUALITY group). + * \param[in] config - Definition of the particular problem. + */ + void SetVolumeOutputFieldsFlameMeshQuality(const CConfig* config); + /*! * \brief Set all scalar (turbulence/species) volume field values for a point. * \param[in] config - Definition of the particular problem. diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index 29032b003026..0f94f033a0ef 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/SU2_CFD/include/output/COutput.hpp @@ -98,6 +98,15 @@ class COutput { string historyFilename; /*!< \brief The history filename*/ ofstream histFile; /*!< \brief Output file stream for the history */ + /*! \brief A dedicated output file for all "Probe"-type custom outputs that share one [x,y,z] location. */ + struct ProbeHistoryFile { + std::vector coords; /*!< \brief Target coordinates (as parsed), used to group outputs at the same point. */ + std::vector names; /*!< \brief Names of all custom outputs probed at this location. */ + ofstream file; /*!< \brief Output file stream, one row appended per history write. */ + }; + /*! \brief One file per unique probe location in CUSTOM_OUTPUTS (MASTER_NODE only, populated in PrepareProbeHistoryFiles). */ + std::vector probeHistoryFiles; + bool cauchyTimeConverged; /*! \brief: Flag indicating that solver is already converged. Needed for writing restart files. */ bool maxTimeDelayActive; /*! \brief: Flag for delaying stop at max_time with 2nd order time stepping. */ @@ -837,6 +846,20 @@ class COutput { */ void PrepareHistoryFile(CConfig *config); + /*! + * \brief Open one dedicated file per "Probe"-type CUSTOM_OUTPUTS entry and write its header. + * MASTER_NODE only, called once during history output preprocessing. + * \param[in] config - Definition of the particular problem. + */ + void PrepareProbeHistoryFiles(const CConfig *config); + + /*! + * \brief Append one row (iteration indices + value) to each dedicated probe history file. + * MASTER_NODE only. Reads already-computed history field values, no extra solver work. + * \param[in] config - Definition of the particular problem. + */ + void SetProbeHistoryFileOutput(const CConfig *config); + /*! * \brief Load up the values of the requested volume fields into ::Local_Data array. * \param[in] config - Definition of the particular problem. @@ -965,6 +988,17 @@ class COutput { * \param[in] solver - The container holding all solution data. * \param[in] iPoint - Index of the point. */ + /*! + * \brief Prepare per-write quantities before the volume data point loop runs. + * Anything global (a reduction over ranks, a sweep over all points) belongs here rather than in + * LoadVolumeData: that is called once per point, so a rank owning no points would never execute + * it, and a collective placed inside it would hang every other rank. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - The container holding all solution data. + */ + inline virtual void PrepareVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver){} + inline virtual void LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){} /*! diff --git a/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp b/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp index fda8056fcd55..0dca2f0f4a95 100644 --- a/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp +++ b/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp @@ -37,8 +37,18 @@ class CSpeciesFlameletVariable final : public CSpeciesVariable { protected: 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 */ + MatrixType source_pd; /*!< \brief PD closure source terms S_{phi_k} (Eq. 16) per control variable, for visualization. */ + MatrixType pd_flux_coeff; /*!< \brief Eq. (14) preferential diffusion flux coefficients (SOURCE_TERM method): per + control variable k, columns [k*pd_terms_per_cv + s] hold the molecular coefficients + D_{phi_k,s} for major species s, and column [k*pd_terms_per_cv + n_major_species] + holds the thermal (Soret) coefficient D^T_{phi_k}. */ + unsigned short pd_terms_per_cv = 0; /*!< \brief Stride of pd_flux_coeff: one molecular coefficient per configured + major species plus the single thermal coefficient. */ su2vector 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. */ + su2activevector hull_miss_dcv1_; /*!< \brief Signed CV1 deviation (query minus nearest hull node) at worst-miss Z level. */ + su2activevector hull_miss_dcv2_; /*!< \brief Signed CV2 deviation (query minus nearest hull node) at worst-miss Z level. */ + su2activevector z_level_dist_; /*!< \brief Distance to nearest table Z level in physical Z units. */ public: /*! @@ -85,10 +95,61 @@ class CSpeciesFlameletVariable final : public CSpeciesVariable { */ inline const su2double* GetScalarLookups(unsigned long iPoint) const override { return lookup_scalar[iPoint]; } + /*! + * \brief Store the PD closure source term S_{phi_k} for control variable iCV (for visualization). + */ + inline void SetScalarSourcePD(unsigned long iPoint, unsigned short iCV, su2double val) { + source_pd(iPoint, iCV) = val; + } + + /*! + * \brief Get the PD closure source terms S_{phi_k} for all control variables at iPoint. + */ + inline const su2double* GetScalarSourcesPD(unsigned long iPoint) const override { return source_pd[iPoint]; } + + /*! + * \brief Store an Eq. (14) preferential diffusion flux coefficient (SOURCE_TERM method). + * \param[in] iPoint - Node index. + * \param[in] iCV - Control variable index. + * \param[in] iTerm - Major species index (molecular), or n_major_species for the thermal coefficient. + * \param[in] val - Coefficient value from the manifold. + */ + inline void SetPDFluxCoeff(unsigned long iPoint, unsigned short iCV, unsigned short iTerm, su2double val) { + pd_flux_coeff(iPoint, iCV * pd_terms_per_cv + iTerm) = val; + } + + /*! + * \brief Get an Eq. (14) preferential diffusion flux coefficient (SOURCE_TERM method). + * \param[in] iPoint - Node index. + * \param[in] iCV - Control variable index. + * \param[in] iTerm - Major species index (molecular), or n_major_species for the thermal coefficient. + */ + inline su2double GetPDFluxCoeff(unsigned long iPoint, unsigned short iCV, unsigned short iTerm) const override { + return pd_flux_coeff(iPoint, iCV * pd_terms_per_cv + iTerm); + } + + /*! + * \brief The whole Eq. (14) coefficient matrix, for the edge flux to gather from. Column + * iCV * GetPDTermsPerCV() + iTerm, with iTerm == n_major_species the thermal coefficient. + */ + inline const MatrixType& GetPDFluxCoeffs() const { return pd_flux_coeff; } + + /*! + * \brief Column stride of GetPDFluxCoeffs(): one coefficient per major species plus the thermal one. + */ + inline unsigned short GetPDTermsPerCV() const { return pd_terms_per_cv; } + 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]; } + inline void SetHullMissDevCV1(unsigned long iPoint, su2double val) override { hull_miss_dcv1_(iPoint) = val; } + inline su2double GetHullMissDevCV1(unsigned long iPoint) const override { return hull_miss_dcv1_(iPoint); } + inline void SetHullMissDevCV2(unsigned long iPoint, su2double val) override { hull_miss_dcv2_(iPoint) = val; } + inline su2double GetHullMissDevCV2(unsigned long iPoint) const override { return hull_miss_dcv2_(iPoint); } + inline void SetZLevelDist(unsigned long iPoint, su2double val) override { z_level_dist_(iPoint) = val; } + inline su2double GetZLevelDist(unsigned long iPoint) const override { return z_level_dist_(iPoint); } + /*! * \brief Store the consumption-rate Jacobian dS_aux_i/dY_aux_i = source_cons_i for user scalar i_aux. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 739a1af75fa5..c589eee4a4fa 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2441,4 +2441,22 @@ class CVariable { inline virtual su2double GetHbyACorrection(unsigned long iPoint, unsigned short iDim) { return 0.0; } inline virtual void SetHbyACorrection(unsigned long iPoint, unsigned short iDim, su2double val_HbyAcorrection) { } + inline virtual const su2double *GetScalarSourcesPD(unsigned long iPoint) const { return nullptr; } + + /*! + * \brief Get an Eq. (14) preferential diffusion flux coefficient (flamelet SOURCE_TERM method). + * \param[in] iPoint - Node index. + * \param[in] iCV - Control variable index. + * \param[in] iTerm - Major species index (molecular), or n_major_species for the thermal coefficient. + */ + inline virtual su2double GetPDFluxCoeff(unsigned long iPoint, unsigned short iCV, unsigned short iTerm) const { + return 0.0; + } + + inline virtual void SetHullMissDevCV1(unsigned long iPoint, su2double val) {} + inline virtual su2double GetHullMissDevCV1(unsigned long iPoint) const { return 0.0; } + inline virtual void SetHullMissDevCV2(unsigned long iPoint, su2double val) {} + inline virtual su2double GetHullMissDevCV2(unsigned long iPoint) const { return 0.0; } + inline virtual void SetZLevelDist(unsigned long iPoint, su2double val) {} + inline virtual su2double GetZLevelDist(unsigned long iPoint) const { return 0.0; } }; diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index dc9ef2ba2b2f..361c892a612f 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -103,7 +103,15 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati PreprocessLookUp(config); if (rank == MASTER_NODE) { - cout << "Preferential diffusion: " << (preferential_diffusion ? "Enabled" : "Disabled") << endl; + if (preferential_diffusion) { + const auto method = flamelet_options.pd_method; + const char* method_str = (method == FLAMELET_PD_METHOD::BETA_CORRECTION) ? "BETA_CORRECTION" + : (method == FLAMELET_PD_METHOD::SOURCE_TERM) ? "SOURCE_TERM" + : "COMBINED"; + cout << "Preferential diffusion: Enabled (" << method_str << ")" << endl; + } else { + cout << "Preferential diffusion: Disabled" << endl; + } } } @@ -212,19 +220,61 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { varnames_LookUp[iLookup] = flamelet_options.lookup_names[iLookup]; } - /*--- Preferential diffusion scalars ---*/ - varnames_PD.resize(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); - val_vars_PD.resize(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); + /*--- Preferential diffusion scalars: only the terms required by the active PD method are looked up. + * + * BETA_CORRECTION: [Beta_ProgVar, Beta_Enth_Thermal, Beta_Enth, Beta_MixFrac] + * + * SOURCE_TERM (model B2 of Schepers & van Oijen, C&F 280 (2025) 114332), layout: + * [0, n_CV) : S_PV, S_h, (S_Z) — Eq. (16) closure sources of the + * NON-major species only, volumetric (rho-weighted) units: + * [kg/(m^3 s)] for PV and Z, [W/m^3] for h. + * [n_CV, n_CV+3) : Y-H2, Y-H2O, Y-H — major species mass fractions [-]. + * [n_CV+3, n_CV+3+3*n_CV) : D_{cv}_{sp} — molecular flux coefficients + * D_{phi_k,i} = c_{phi_k,i} (rho*D_i - lambda/cp) of Eq. (14), + * cv in {PV, h, Z} (CV order), sp in {H2, H2O, H}; + * units [kg/(m s)] for PV/Z rows, [W/m] for the h row + * (c_{h,i} = h_i(T) is folded in at generation time). + * [n_CV+3+3*n_CV, ...) : DT_PV, DT_h, (DT_Z) — combined Soret coefficients + * DT_{phi_k} = sum_sp c_{phi_k,sp} DT_sp, with DT_sp the + * Cantera thermal-diffusion coefficient [kg/(m s)] (species + * flux j_sp = -(rho D_sp grad Y_sp + DT_sp/T grad T)); units + * [kg/(m s)] for PV/Z, [W/m] for h. The 1/T factor is NOT + * tabulated; it is applied at runtime with the local CFD T. + * The runtime flux applied in CSpeciesFlameletSolver::Viscous_Residual is + * J_{phi_k} = - sum_sp D_{cv}_{sp} grad(Y_sp) - (DT_{cv}/T) grad(T), + * entering the transport equation as -div(J). ---*/ + const auto pd_method = flamelet_options.pd_method; + const bool use_beta = (pd_method != FLAMELET_PD_METHOD::SOURCE_TERM); + const bool use_src = (pd_method != FLAMELET_PD_METHOD::BETA_CORRECTION); + const unsigned n_beta_vars = use_beta ? FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS : 0u; + n_pd_major_species = use_src ? flamelet_options.n_pd_major_species : 0u; + const unsigned n_majors = n_pd_major_species; + const unsigned n_src_vars = use_src ? (n_control_vars + n_majors + (n_majors + 1) * n_control_vars) : 0u; + const unsigned n_pd_terms = n_beta_vars + n_src_vars; + varnames_PD.resize(n_pd_terms); + val_vars_PD.resize(n_pd_terms, 0.0); + + if (use_beta) { + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR] = "Beta_ProgVar"; + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL] = "Beta_Enth_Thermal"; + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH] = "Beta_Enth"; + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC] = "Beta_MixFrac"; + } - varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR] = "Beta_ProgVar"; - varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL] = "Beta_Enth_Thermal"; - varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH] = "Beta_Enth"; - varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC] = "Beta_MixFrac"; + if (use_src) { + const auto* major_names = flamelet_options.pd_major_species_names; + const auto* cv_names = flamelet_options.controlling_variable_names; - val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR] = beta_progvar; - val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL] = beta_enth_thermal; - val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH] = beta_enth; - val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC] = beta_mixfrac; + for (auto iCV = 0u; iCV < n_control_vars; iCV++) + varnames_PD[n_beta_vars + iCV] = "Res_" + cv_names[iCV]; + + unsigned idx = n_beta_vars + n_control_vars; + for (auto iSp = 0u; iSp < n_majors; iSp++) varnames_PD[idx++] = "Y-" + major_names[iSp]; + for (auto iCV = 0u; iCV < n_control_vars; iCV++) + for (auto iSp = 0u; iSp < n_majors; iSp++) + varnames_PD[idx++] = "D_" + cv_names[iCV] + "_" + major_names[iSp]; + for (auto iCV = 0u; iCV < n_control_vars; iCV++) varnames_PD[idx++] = "DT_" + cv_names[iCV]; + } preferential_diffusion = flamelet_options.preferential_diffusion; @@ -267,7 +317,18 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { LUT_idx_LookUp.push_back(LUT_idx); } if (preferential_diffusion) { - for (auto iVar = 0u; iVar < varnames_PD.size(); iVar++) { + /*--- Check all preferential diffusion variables in one pass, so that a manifold missing + several of the required columns (e.g. the Eq. (14) flux coefficients of the SOURCE_TERM + method) reports the complete list instead of aborting on the first one. ---*/ + std::string missing_vars; + for (const auto& name : varnames_PD) + if (!look_up_table->CheckForVariables({name})) missing_vars += "\n " + name; + if (!missing_vars.empty()) + SU2_MPI::Error("The manifold is missing the following variables required by the active " + "PREFERENTIAL_DIFFUSION_METHOD (see CFluidFlamelet::PreprocessLookUp for " + "their definitions and units):" + missing_vars, CURRENT_FUNCTION); + + for (auto iVar=0u; iVar < varnames_PD.size(); iVar++) { LUT_idx_PD.push_back(look_up_table->GetIndexOfVar(varnames_PD[iVar])); } } @@ -342,3 +403,44 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca AD::EndPreacc(); return extrapolation; } + +void CFluidFlamelet::GetTableCVBounds(su2double& cv1_min, su2double& cv1_max, + su2double& cv2_min, su2double& cv2_max, + su2double& cv3_min, su2double& cv3_max) const { + if (look_up_table == nullptr) { + cv1_min = cv2_min = cv3_min = 0.0; + cv1_max = cv2_max = cv3_max = 1.0; + return; + } + auto lx0 = look_up_table->GetTableLimitsX(0); + cv1_min = *lx0.first; cv1_max = *lx0.second; + auto ly0 = look_up_table->GetTableLimitsY(0); + cv2_min = *ly0.first; cv2_max = *ly0.second; + for (unsigned long l = 1; l < look_up_table->GetNTableLevels(); ++l) { + auto lx = look_up_table->GetTableLimitsX(l); + if (*lx.first < cv1_min) cv1_min = *lx.first; + if (*lx.second > cv1_max) cv1_max = *lx.second; + auto ly = look_up_table->GetTableLimitsY(l); + if (*ly.first < cv2_min) cv2_min = *ly.first; + if (*ly.second > cv2_max) cv2_max = *ly.second; + } + auto lz = look_up_table->GetTableLimitsZ(); + cv3_min = lz.first; + cv3_max = lz.second; +} + +void CFluidFlamelet::ResetHullMissDistance() { + if (look_up_table) look_up_table->ResetHullMissDistance(); +} + +su2double CFluidFlamelet::GetHullMissCV1Dev() const { + return look_up_table ? look_up_table->GetHullMissCV1Dev() : 0.0; +} + +su2double CFluidFlamelet::GetHullMissCV2Dev() const { + return look_up_table ? look_up_table->GetHullMissCV2Dev() : 0.0; +} + +su2double CFluidFlamelet::GetDistanceToNearestZLevel(su2double val_CV3) const { + return look_up_table ? look_up_table->GetDistanceToNearestZLevel(val_CV3) : 0.0; +} diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 2ac2693b9952..ae8e8e07ac78 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -755,7 +755,8 @@ void CFlowOutput::ConvertVariableSymbolsToIndices(const CPrimitiveIndices& nameToIndex, const std::string& var) { /*--- Primitives of the flow solver. ---*/ @@ -769,6 +770,18 @@ void CFlowOutput::ConvertVariableSymbolsToIndices(const CPrimitiveIndicesGetPrimitive(iPoint, varIdx); } + if (solIdx == SPECIES_SOL && varIdx >= 8) { + /*--- Sub-ranges within the species solver's slot select scalar sources / PD closure + sources / PD beta scalars instead of the solution (see ConvertVariableSymbolsToIndices). + All of these are already computed earlier this iteration in the species solver's + Preprocessing, so this is just an array read, no extra solver work. ---*/ + if (varIdx >= 24) { + return solver[SPECIES_SOL]->GetNodes()->GetAuxVar(iPoint, varIdx - 24); + } + if (varIdx >= 16) { + const auto* src_pd = solver[SPECIES_SOL]->GetNodes()->GetScalarSourcesPD(iPoint); + return src_pd ? src_pd[varIdx - 16] : su2double(0.0); + } + const auto* src = solver[SPECIES_SOL]->GetNodes()->GetScalarSources(iPoint); + return src ? src[varIdx - 8] : su2double(0.0); + } return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx); } return *output.otherOutputs[i - CustomOutput::NOT_A_VARIABLE]; @@ -1567,6 +1595,14 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSource(const CConfig* config) { if (cv_source_name.compare("NULL") != 0) AddVolumeOutput("SOURCE_"+cv_name, "Source_" + cv_name, "SOURCE", "Source " + cv_name); } + /*--- PD closure source terms S_{phi_k} from Eq. 16 (Schepers & van Oijen 2025), for debugging. ---*/ + if (flamelet_config_options.preferential_diffusion && + flamelet_config_options.pd_method == FLAMELET_PD_METHOD::SOURCE_TERM) { + for (auto iCV = 0u; iCV < flamelet_config_options.n_control_vars; iCV++) { + const auto& cv_name = flamelet_config_options.controlling_variable_names[iCV]; + AddVolumeOutput("PD_SOURCE_" + cv_name, "PD_Source_" + cv_name, "SOURCE", "PD closure source S_{" + cv_name + "} (Eq. 16)"); + } + } /*--- no source term for enthalpy ---*/ /*--- auxiliary species source terms ---*/ for (auto iReactant=0u; iReactantGetKind_Species_Model() != SPECIES_MODEL::FLAMELET) return; + + const auto* Node_Species = solver[SPECIES_SOL]->GetNodes(); + su2double pv_loc_min = std::numeric_limits::max(); + su2double pv_loc_max = std::numeric_limits::lowest(); + for (auto iPt = 0ul; iPt < geometry->GetnPointDomain(); ++iPt) { + const su2double pv = Node_Species->GetSolution(iPt, 0); + pv_loc_min = std::min(pv_loc_min, pv); + pv_loc_max = std::max(pv_loc_max, pv); + } + su2double pv_global_min, pv_global_max; + SU2_MPI::Allreduce(&pv_loc_min, &pv_global_min, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&pv_loc_max, &pv_global_max, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); + flamelet_pv_range = max(pv_global_max - pv_global_min, 1e-10); +} + +void CFlowOutput::SetVolumeOutputFieldsFlameMeshQuality(const CConfig* config) { + if (config->GetKind_Species_Model() != SPECIES_MODEL::FLAMELET) return; + const auto& fco = config->GetFlameletParsedOptions(); + + AddVolumeOutput("FLAME_GRAD_C", "grad_C_mag", "FLAME_QUALITY", "|grad(C)| progress variable gradient magnitude [1/m]"); + AddVolumeOutput("FLAME_C_PLUS", "C_Plus", "FLAME_QUALITY", "Flame resolution index C+ = |grad(C)|*h_cell/dC_range; target < 0.067 (>=15 cells), flag > 0.2"); + AddVolumeOutput("FLAME_CHI_C", "chi_C", "FLAME_QUALITY", "Scalar dissipation rate chi_C = 2*D_C*|grad(C)|^2 [1/s]"); + AddVolumeOutput("FLAME_DA_LOCAL", "Da_local", "FLAME_QUALITY", "Local Damkohler number Da = h_cell*omega_C / (|u|*C)"); + if (fco.preferential_diffusion && fco.pd_method != FLAMELET_PD_METHOD::SOURCE_TERM) + AddVolumeOutput("FLAME_GRAD_BETA", "grad_beta_C_mag", "FLAME_QUALITY", "|grad(beta_C)| Schwab-Zeldovich correction scalar gradient [1/m]"); } void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* const* solver, const CGeometry* geometry, @@ -1757,6 +1833,15 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con if (source_name.compare("NULL") != 0) SetVolumeOutputValue("SOURCE_" + cv_name, iPoint, Node_Species->GetScalarSources(iPoint)[iCV]); } + /*--- PD closure source terms S_{phi_k} (Eq. 16) for debugging. ---*/ + if (flamelet_config_options.preferential_diffusion && + flamelet_config_options.pd_method == FLAMELET_PD_METHOD::SOURCE_TERM) { + const auto* pd_src = Node_Species->GetScalarSourcesPD(iPoint); + for (auto iCV = 0u; iCV < flamelet_config_options.n_control_vars; iCV++) { + const auto& cv_name = flamelet_config_options.controlling_variable_names[iCV]; + SetVolumeOutputValue("PD_SOURCE_" + cv_name, iPoint, pd_src[iCV]); + } + } /*--- auxiliary species transport equations ---*/ for (unsigned short i_scalar=0; i_scalarGetScalarLookups(iPoint)[i_lookup]); } - SetVolumeOutputValue("TABLE_MISSES", iPoint, Node_Species->GetTableMisses(iPoint)); + SetVolumeOutputValue("TABLE_MISSES", iPoint, Node_Species->GetTableMisses(iPoint)); + SetVolumeOutputValue("HULL_MISS_DELTA_CV1", iPoint, Node_Species->GetHullMissDevCV1(iPoint)); + SetVolumeOutputValue("HULL_MISS_DELTA_CV2", iPoint, Node_Species->GetHullMissDevCV2(iPoint)); + SetVolumeOutputValue("Z_LEVEL_DIST", iPoint, Node_Species->GetZLevelDist(iPoint)); + + /*--- Flame mesh quality diagnostics (FLAME_QUALITY group). The global progress variable + range they normalise against is computed once per write in PrepareVolumeData. ---*/ + { + /*--- Cell characteristic length h_cell = V^(1/nDim). ---*/ + const su2double vol = geometry->nodes->GetVolume(iPoint); + const su2double h_cell = pow(vol, 1.0 / nDim); + + /*--- |grad(C)| ---*/ + su2double grad_C_sq = 0.0; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) { + const su2double gc = Node_Species->GetGradient(iPoint, 0, iDim); + grad_C_sq += gc * gc; + } + const su2double grad_C_mag = sqrt(grad_C_sq); + + SetVolumeOutputValue("FLAME_GRAD_C", iPoint, grad_C_mag); + SetVolumeOutputValue("FLAME_C_PLUS", iPoint, grad_C_mag * h_cell / flamelet_pv_range); + SetVolumeOutputValue("FLAME_CHI_C", iPoint, 2.0 * Node_Species->GetDiffusivity(iPoint, 0) * grad_C_sq); + + /*--- Local Damköhler: Da = h_cell * omega_C / (|u| * C) ---*/ + const su2double C_pv = Node_Species->GetSolution(iPoint, 0); + const su2double omega_C = Node_Species->GetScalarSources(iPoint)[0]; + su2double vel_sq = 0.0; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) { + const su2double ui = Node_Flow->GetVelocity(iPoint, iDim); + vel_sq += ui * ui; + } + SetVolumeOutputValue("FLAME_DA_LOCAL", iPoint, (h_cell * omega_C) / max(sqrt(vel_sq) * C_pv, 1e-10)); + + /*--- |grad(beta_C)| for PD-BETA_CORRECTION only ---*/ + if (flamelet_config_options.preferential_diffusion && + flamelet_config_options.pd_method != FLAMELET_PD_METHOD::SOURCE_TERM) { + su2double grad_beta_sq = 0.0; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) { + const su2double gb = Node_Species->GetAuxVarGradient(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR, iDim); + grad_beta_sq += gb * gb; + } + SetVolumeOutputValue("FLAME_GRAD_BETA", iPoint, sqrt(grad_beta_sq)); + } + } } break; diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp index c9ca947546cb..919359fdd8b0 100644 --- a/SU2_CFD/src/output/COutput.cpp +++ b/SU2_CFD/src/output/COutput.cpp @@ -26,6 +26,8 @@ */ #include +#include +#include #include #include "../../../Common/include/geometry/CGeometry.hpp" @@ -289,7 +291,10 @@ void COutput::OutputScreenAndHistory(CConfig *config) { if (rank == MASTER_NODE && !noWriting) { - if (WriteHistoryFileOutput(config)) SetHistoryFileOutput(config); + if (WriteHistoryFileOutput(config)) { + SetHistoryFileOutput(config); + SetProbeHistoryFileOutput(config); + } if (WriteScreenHeader(config)) SetScreenHeader(config); @@ -1287,8 +1292,10 @@ void COutput::PreprocessHistoryOutput(CConfig *config, bool wrt){ if (rank == MASTER_NODE && !noWriting){ /*--- Open history file and print the header ---*/ - if (!config->GetMultizone_Problem() || config->GetWrt_ZoneHist()) + if (!config->GetMultizone_Problem() || config->GetWrt_ZoneHist()) { PrepareHistoryFile(config); + PrepareProbeHistoryFiles(config); + } total_width = nRequestedScreenFields*fieldWidth + (nRequestedScreenFields-1); @@ -1396,6 +1403,63 @@ void COutput::PrepareHistoryFile(CConfig *config){ } +void COutput::PrepareProbeHistoryFiles(const CConfig* config) { + + /*--- Derive a common base name/extension from the main history file so probe files + land alongside it and follow the same zone/restart-iteration naming. ---*/ + + string ext = ".csv"; + if (config->GetTabular_FileFormat() == TAB_OUTPUT::TAB_TECPLOT) ext = ".dat"; + + string base = historyFilename; + PrintingToolbox::TrimExtension(ext, base); + + probeHistoryFiles.clear(); + + /*--- Group "Probe"-type custom outputs that target the same [x,y,z] location into a single + dedicated file, so every variable probed at one physical point lands together (e.g. + Probe{SPECIES[0]}, Probe{SPECIES[2]}, Probe{VELOCITY_Y} all at the same coordinates). ---*/ + for (const auto& output : customOutputs) { + if (output.type != OperationType::PROBE) continue; + + auto it = std::find_if(probeHistoryFiles.begin(), probeHistoryFiles.end(), + [&](const ProbeHistoryFile& f) { return f.coords == output.markers; }); + if (it == probeHistoryFiles.end()) { + probeHistoryFiles.emplace_back(); + it = std::prev(probeHistoryFiles.end()); + it->coords = output.markers; + } + it->names.push_back(output.name); + } + + for (auto& probeFile : probeHistoryFiles) { + string suffix; + for (const auto& name : probeFile.names) suffix += "_" + name; + probeFile.file.open(base + "_probe" + suffix + ext, ios::out); + + probeFile.file << "\"Time_Iter\",\"Outer_Iter\",\"Inner_Iter\""; + for (const auto& name : probeFile.names) probeFile.file << ",\"" << name << "\""; + probeFile.file << endl; + } +} + +void COutput::SetProbeHistoryFileOutput(const CConfig* config) { + + /*--- Probe values were already computed (and Allreduced) earlier this iteration as part of the + regular history output evaluation (CFlowOutput::SetCustomOutputs). This only re-reads the + cached value(s) from the history map, so the extra cost per probe location is one map + lookup per variable. ---*/ + + for (auto& probeFile : probeHistoryFiles) { + probeFile.file << curTimeIter << ", " << curOuterIter << ", " << curInnerIter; + for (const auto& name : probeFile.names) { + probeFile.file << ", " << std::setprecision(config->GetOutput_Precision()) << GetHistoryFieldValue(name); + } + probeFile.file << "\n"; + probeFile.file.flush(); + } +} + void COutput::CheckHistoryOutput(unsigned short nZone) { /*--- Set screen convergence output header and remove unavailable fields ---*/ @@ -1735,6 +1799,10 @@ void COutput::LoadDataIntoSorter(CConfig* config, CGeometry* geometry, CSolver** curGetFieldIndex = 0; fieldGetIndexCache.clear(); + /*--- Per-write preparation. Runs unconditionally on every rank, which is what makes it safe for + collectives, unlike anything reached from the per-point loop below. ---*/ + PrepareVolumeData(config, geometry, solver); + if (femOutput) { /*--- Create an object of the class CMeshFEM_DG and retrieve the necessary diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 1d39f9263381..bc39cd6e44f6 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -81,7 +81,8 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { SU2_ZONE_SCOPED - unsigned long n_not_in_domain_global = 0; + unsigned long n_not_in_domain_local = 0, n_not_in_domain_global = 0; + unsigned long n_miss_pv = 0, n_miss_enth = 0, n_miss_mf = 0, n_miss_hull = 0; vector scalars_vector(nVar); unsigned long spark_iter_start, spark_duration; @@ -130,6 +131,14 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver /* Flame thickness correction factors */ su2double F{1.0}, F_source{1.0}; + /*--- Pre-fetch CV bounds for verbose miss classification (read-only, safe for all threads). ---*/ + const bool verbose_misses = flamelet_config_options.verbose_misses; + const bool has_mf = (flamelet_config_options.n_control_vars == 3); + su2double cv1_min = 0, cv1_max = 1, cv2_min = 0, cv2_max = 1, cv3_min = 0, cv3_max = 1; + if (verbose_misses) + static_cast(solver_container[FLOW_SOL]->GetFluidModel()) + ->GetTableCVBounds(cv1_min, cv1_max, cv2_min, cv2_max, cv3_min, cv3_max); + SU2_OMP_FOR_STAT(omp_chunk_size) for (auto i_point = 0u; i_point < nPoint; i_point++) { CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); @@ -143,6 +152,9 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar]; + /*--- Reset hull-miss distance accumulator for this point before any table lookups. ---*/ + fluid_model_local->ResetHullMissDistance(); + /*--- Only apply thickened flame correction factor to sources for steady problems. ---*/ unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector, F_source); @@ -164,6 +176,19 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver nodes->SetTableMisses(i_point, misses); SU2_OMP_ATOMIC n_not_in_domain_local += misses; + if (verbose_misses && misses != 0) { + const su2double C = scalars_vector[FLAMELET_SCALAR_VARIABLES::I_PROGVAR]; + const su2double h = scalars_vector[FLAMELET_SCALAR_VARIABLES::I_ENTH]; + if (has_mf && (scalars_vector[FLAMELET_SCALAR_VARIABLES::I_MIXFRAC] < cv3_min || + scalars_vector[FLAMELET_SCALAR_VARIABLES::I_MIXFRAC] > cv3_max)) + n_miss_mf++; + else if (C < cv1_min || C > cv1_max) + n_miss_pv++; + else if (h < cv2_min || h > cv2_max) + n_miss_enth++; + else + n_miss_hull++; + } /*--- Obtain passive look-up scalars. ---*/ SetScalarLookUps(fluid_model_local, i_point, scalars_vector); @@ -179,6 +204,17 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver if (flamelet_config_options.preferential_diffusion) SetPreferentialDiffusionScalars(fluid_model_local, i_point, scalars_vector); + /*--- Store signed CV deviations (query minus nearest hull node) at the worst-miss Z level. ---*/ + nodes->SetHullMissDevCV1(i_point, fluid_model_local->GetHullMissCV1Dev()); + nodes->SetHullMissDevCV2(i_point, fluid_model_local->GetHullMissCV2Dev()); + + /*--- Store distance to nearest table Z level (zero for 2D tables). ---*/ + if (has_mf) { + const su2double z_dist = static_cast(fluid_model_local) + ->GetDistanceToNearestZLevel(scalars_vector[I_MIXFRAC]); + nodes->SetZLevelDist(i_point, z_dist); + } + if (!Output) LinSysRes.SetBlock_Zero(i_point); } END_SU2_OMP_FOR @@ -189,19 +225,62 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver if ((rank == MASTER_NODE) && (n_not_in_domain_global > 0)) cout << "Number of points outside manifold domain: " << n_not_in_domain_global << endl;) - /*--- Compute preferential diffusion scalar gradients. ---*/ - if (flamelet_config_options.preferential_diffusion) { - switch (config->GetKind_Gradient_Method()) { - case GREEN_GAUSS: - SetAuxVar_Gradient_GG(geometry, config); - break; - case WEIGHTED_LEAST_SQUARES: - SetAuxVar_Gradient_LS(geometry, config); - break; - default: - break; + if (verbose_misses) { + unsigned long miss_cats[4] = {n_miss_pv, n_miss_enth, n_miss_mf, n_miss_hull}; + unsigned long miss_cats_global[4] = {}; + SU2_MPI::Reduce(miss_cats, miss_cats_global, 4, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm()); + if ((rank == MASTER_NODE) && (n_not_in_domain_global > 0)) + cout << " Progress variable: " << miss_cats_global[0] + << " | Total Enthalpy: " << miss_cats_global[1] + << " | Mixture Fraction: " << miss_cats_global[2] + << " | Hull: " << miss_cats_global[3] << endl; + } + + /*--- Compute auxiliary-variable gradients: β scalars for the BETA_CORRECTION method, the major + species mass fractions (PREFERENTIAL_DIFFUSION_MAJOR_SPECIES) for the SOURCE_TERM method. These + gradients ARE the multi-dimensional grad(Y_i) of Eq. (14) — the whole point of tabulating the + coefficients per species rather than pre-contracting them against the 1D flamelet gradients — so + they must exist for every gradient method. LEAST_SQUARES used to fall through the switch, leaving + the auxiliary gradients at zero and silently disabling preferential diffusion altogether; it is + routed to the least-squares path instead. The companion grad(T) of Eq. (14) is the flow solver's + primitive temperature gradient, which is always available for a viscous run. ---*/ + /*--- Zero the Eq. (14) thermal (Soret) coefficients on viscous walls, matching the reference + implementation. The justification is manifold support, not impermeability: the tabulated Soret + coefficients come from freely propagating and burner-stabilized flamelets, none of which has a + quench layer, so applying them against the steepest temperature gradient in the domain uses the + manifold outside the states it was built from. Note that this does NOT remove a flux through the + wall face -- Viscous_Residual visits interior edges only, so the wall face never carries an + Eq. (14) flux. What it changes is the near-wall INTERIOR edges, because the scalar numerics + averages the two nodal coefficients. The molecular coefficients are deliberately left alone: + they act on the major species gradients, which the wall boundary condition already controls. + This is a modelling choice with a measurable consequence -- on the burner case it moves the wall + heat flux by up to 1.2% and the wall temperature by 0.4 K, decaying to nothing within 0.5 mm. ---*/ + if (flamelet_config_options.preferential_diffusion && + flamelet_config_options.pd_method == FLAMELET_PD_METHOD::SOURCE_TERM) { + auto* wall_nodes = static_cast(nodes); + const auto n_CV_wall = flamelet_config_options.n_control_vars; + const auto i_thermal = FlameletPDThermalTerm(flamelet_config_options.n_pd_major_species); + + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { + /*--- Viscous walls only. GetSolid_Wall would also match EULER_WALL, a slip boundary with no + quench layer and so nothing to suppress. ---*/ + if (!config->GetViscous_Wall(iMarker)) continue; + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) { + const auto iPoint_wall = geometry->vertex[iMarker][iVertex]->GetNode(); + for (auto iCV = 0u; iCV < n_CV_wall; iCV++) + wall_nodes->SetPDFluxCoeff(iPoint_wall, iCV, i_thermal, 0.0); + } + END_SU2_OMP_FOR } } + + if (flamelet_config_options.preferential_diffusion) { + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) + SetAuxVar_Gradient_GG(geometry, config); + else + SetAuxVar_Gradient_LS(geometry, config); + } /*--- Clear Residual and Jacobian. Upwind second order reconstruction and gradients ---*/ CommonPreprocessing(geometry, config, Output); } @@ -709,14 +788,87 @@ unsigned long CSpeciesFlameletSolver::SetPreferentialDiffusionScalars(CFluidMode unsigned long iPoint, const vector& scalars) { SU2_ZONE_SCOPED - /*--- Retrieve the preferential diffusion scalar values from the manifold. ---*/ + /*--- Retrieve preferential diffusion scalars from the manifold. + * BETA_CORRECTION: [Beta_ProgVar, Beta_Enth_Thermal, Beta_Enth, Beta_MixFrac] + * SOURCE_TERM: [Res_] + [Y-] + Eq. (14) flux coefficients; + * the full layout is documented in CFluidFlamelet::PreprocessLookUp and must stay in + * sync with the sizing below. ---*/ + const auto pd_method = flamelet_config_options.pd_method; + const bool use_beta = (pd_method != FLAMELET_PD_METHOD::SOURCE_TERM); + const bool use_src = (pd_method != FLAMELET_PD_METHOD::BETA_CORRECTION); + const auto n_CV = flamelet_config_options.n_control_vars; + const unsigned n_majors = flamelet_config_options.n_pd_major_species; + const unsigned n_beta_vars = use_beta ? FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS : 0u; + const unsigned n_pd_terms = n_beta_vars + (use_src ? (n_CV + n_majors + (n_majors + 1) * n_CV) : 0u); + vector pref_diff_scalar(n_pd_terms); + unsigned long misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::PREFDIF, pref_diff_scalar); + + /*--- Store beta scalars as aux vars; their gradients are needed for the viscous residual. ---*/ + if (use_beta) { + for (auto i_beta = 0u; i_beta < FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; i_beta++) + nodes->SetAuxVar(iPoint, i_beta, pref_diff_scalar[i_beta]); + } - vector beta_scalar(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); - unsigned long misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::PREFDIF, beta_scalar); + /*--- Major species mass fractions (stored as aux vars, gradients needed for the Eq. (14) + fluxes) and the flux coefficients. The mass fractions are stored unconditionally (their + gradients must remain continuous across the hull boundary), and the runtime fluxes vanish + naturally in composition-uniform regions, so no flame-zone gating is required. + + The flux coefficients, however, are only meaningful ON the manifold: outside the hull the + lookup extrapolates (nearest-neighbor), freezing the coefficients at hull-edge values. At + states driven off the manifold — e.g. wall-quenched burnt gas at an isothermal wall, whose + correct enthalpy h(T_wall) lies below the tabulated h-range — a frozen Soret coefficient + against the quench-layer temperature gradient keeps pumping the state further off-hull + (over-burnt PV, runaway h excursion). Zeroing the coefficients at off-manifold states cuts + this feedback: such nodes fall back to the baseline (unity-Lewis) diffusion, which pulls + them back toward the hull. The deferred-correction stabilization is built from the same + coefficients and vanishes there consistently. ---*/ + if (use_src) { + auto* fn_pd = static_cast(nodes); + unsigned idx = n_beta_vars + n_CV; + for (auto iSp = 0u; iSp < n_majors; iSp++) nodes->SetAuxVar(iPoint, iSp, pref_diff_scalar[idx++]); + const su2double coeff_scale = (misses == 0) ? 1.0 : 0.0; + for (auto iCV = 0u; iCV < n_CV; iCV++) + for (auto iSp = 0u; iSp < n_majors; iSp++) + fn_pd->SetPDFluxCoeff(iPoint, iCV, iSp, coeff_scale * pref_diff_scalar[idx++]); + for (auto iCV = 0u; iCV < n_CV; iCV++) + fn_pd->SetPDFluxCoeff(iPoint, iCV, n_majors, coeff_scale * pref_diff_scalar[idx++]); + } - for (auto i_beta = 0u; i_beta < FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; i_beta++) { - nodes->SetAuxVar(iPoint, i_beta, beta_scalar[i_beta]); + /*--- Accumulate PD closure source terms (Eq. 16, Schepers & van Oijen 2025) into scalar sources. + * Also cache the raw values in source_pd for visualization. + * If the PREFDIF lookup itself misses, the returned values are extrapolated outside the manifold + * and physically meaningless — skip the accumulation to avoid driving the solution further + * off-manifold (see also the THERM lookup misses tracked by SetScalarSources). + * + * Manifold validity is the ONLY gate. An earlier version additionally required chemical + * activity (chem_src_pv > 0), on the reasoning that the manifold carries flamelet-gradient + * structure at every state while the bulk unburnt and burnt regions of the CFD are + * gradient-free. That test did not do what it claimed: SetScalarSources clips the progress + * variable source with fmax(0, .), so the condition only fails where the TABULATED source was + * negative -- a burnt-edge state that occurs inside the flame, not in the bulk. Measured on a + * converged solution it closed at 320 nodes, every one of them in the preheat zone at + * 304-567 K, and stayed open across all 21k gradient-free bulk nodes: the exact inverse of its + * stated purpose, and a one-cell discontinuity in the source field where it did close. The + * preheat zone is also where the preferential diffusion flux divergence is largest for + * hydrogen, so the test suppressed the model precisely where it matters most. ---*/ + if (use_src) { + auto* flamelet_nodes = static_cast(nodes); + if (misses == 0) { + for (auto iCV = 0u; iCV < n_CV; ++iCV) { + const su2double src_pd = pref_diff_scalar[n_beta_vars + iCV]; + flamelet_nodes->SetScalarSourcePD(iPoint, iCV, src_pd); + nodes->SetScalarSource(iPoint, iCV, nodes->GetScalarSources(iPoint)[iCV] + src_pd); + } + } else { + /*--- Off-manifold: no closure source is applied, and the cached value is cleared. Leaving it + alone would let the PD_SOURCE_ visualisation output retain the previous iteration's + value at exactly the nodes where no source was applied. ---*/ + for (auto iCV = 0u; iCV < n_CV; ++iCV) + flamelet_nodes->SetScalarSourcePD(iPoint, iCV, 0.0); + } } + return misses; } @@ -768,6 +920,11 @@ unsigned long CSpeciesFlameletSolver::GetEnthFromTemp(CFluidModel* fluid_model, *val_enth = enth_iter; if (counter >= counter_limit) { + /*--- MASTER_NODE only: this is reached per point, so an unguarded write would be emitted by + every rank for every failing node. ---*/ + if (rank == MASTER_NODE) + cout << " !!! GetEnthFromTemp: Newton iteration did not converge in " << counter_limit + << " iterations, delta_temp_iter = " << delta_temp_iter << endl; exit_code = 1; } @@ -776,17 +933,26 @@ unsigned long CSpeciesFlameletSolver::GetEnthFromTemp(CFluidModel* fluid_model, su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalar_solution, const su2double T_ignition) { SU2_ZONE_SCOPED - su2double scalars[MAXNVAR], delta = 1e-3; + su2double scalars[MAXNVAR], delta = 1e-6; for (auto iVar = 0u; iVar < nVar; iVar++) scalars[iVar] = scalar_solution[iVar]; - bool outside = false; + bool outside = false, hit_t_ignition = false; scalars[I_PROGVAR] += delta; while (!outside) { /*--- 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; + if (fluid_model->GetExtrapolation() == 1) { + outside = true; + } else if (fluid_model->GetTemperature() > T_ignition) { + outside = true; + hit_t_ignition = true; + } scalars[I_PROGVAR] += delta; } - su2double pv_burnt = scalars[I_PROGVAR] - delta; + /*--- When the loop exits via extrapolation the last valid PV is two deltas back + * (the loop increments once after setting outside=true). When it exits via + * the ignition-temperature threshold that state is physically valid and used directly. ---*/ + su2double pv_burnt = scalars[I_PROGVAR] - (hit_t_ignition ? delta : 2 * 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; diff --git a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp index a981655b39af..acd05507be4f 100644 --- a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp @@ -52,11 +52,40 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf, const auto& flamelet_config_options = config->GetFlameletParsedOptions(); source_scalar.resize(nPoint, flamelet_config_options.n_scalars) = su2double(0.0); lookup_scalar.resize(nPoint, flamelet_config_options.n_lookups) = su2double(0.0); + source_pd.resize(nPoint, flamelet_config_options.n_control_vars) = su2double(0.0); table_misses.resize(nPoint) = 0; + hull_miss_dcv1_.resize(nPoint) = su2double(0.0); + hull_miss_dcv2_.resize(nPoint) = su2double(0.0); + z_level_dist_.resize(nPoint) = su2double(0.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); - Grad_AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS, nDim, 0.0); + const bool beta_correction_active = + flamelet_config_options.preferential_diffusion && + flamelet_config_options.pd_method != FLAMELET_PD_METHOD::SOURCE_TERM; + const bool source_term_active = + flamelet_config_options.preferential_diffusion && + flamelet_config_options.pd_method != FLAMELET_PD_METHOD::BETA_CORRECTION; + + /*--- Auxiliary variables: β scalars for the BETA_CORRECTION method, major species mass + fractions (PREFERENTIAL_DIFFUSION_MAJOR_SPECIES) for the SOURCE_TERM method, whose gradients drive the Eq. (14) + preferential diffusion fluxes (Schepers & van Oijen, C&F 280, 2025). Note that nAuxVar must + be set: the SetAuxVar_Gradient_* routines iterate over GetnAuxVar() variables, so leaving it + at its default of zero silently disables the aux-variable gradients. ---*/ + if (beta_correction_active) { + nAuxVar = FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; + } else if (source_term_active) { + nAuxVar = flamelet_config_options.n_pd_major_species; + } + if (nAuxVar > 0) { + AuxVar.resize(nPoint, nAuxVar) = su2double(0.0); + Grad_AuxVar.resize(nPoint, nAuxVar, nDim, 0.0); + } + + /*--- Eq. (14) flux coefficients for the SOURCE_TERM method: per control variable, the + molecular coefficients D_{phi_k,i} for each major species plus one thermal (Soret) + coefficient D^T_{phi_k}. ---*/ + if (source_term_active) { + pd_terms_per_cv = FlameletPDTermsPerCV(flamelet_config_options.n_pd_major_species); + pd_flux_coeff.resize(nPoint, flamelet_config_options.n_control_vars * pd_terms_per_cv) = su2double(0.0); } } diff --git a/SU2_PY/SU2/run/interface.py b/SU2_PY/SU2/run/interface.py index ef1f58f5ca2a..64364475c064 100644 --- a/SU2_PY/SU2/run/interface.py +++ b/SU2_PY/SU2/run/interface.py @@ -30,6 +30,7 @@ # ---------------------------------------------------------------------- import os, sys, shutil, copy +import shlex import subprocess from ..io import Config from ..util import which @@ -46,7 +47,6 @@ ) from exc sys.path.append(SU2_RUN) -quote = '"' if sys.platform == "win32" else "" # SU2 suite run command template base_Command = os.path.join(SU2_RUN, "%s") @@ -98,7 +98,7 @@ def CFD(config): processes = konfig["NUMBER_PART"] - the_Command = "SU2_CFD_DIRECTDIFF%s %s" % (quote, tempname) + the_Command = "SU2_CFD_DIRECTDIFF %s" % tempname elif auto_diff: tempname = "config_CFD_AD.cfg" @@ -106,7 +106,7 @@ def CFD(config): processes = konfig["NUMBER_PART"] - the_Command = "SU2_CFD_AD%s %s" % (quote, tempname) + the_Command = "SU2_CFD_AD %s" % tempname else: tempname = "config_CFD.cfg" @@ -114,7 +114,7 @@ def CFD(config): processes = konfig["NUMBER_PART"] - the_Command = "SU2_CFD%s %s" % (quote, tempname) + the_Command = "SU2_CFD %s" % tempname the_Command = build_command(the_Command, processes) run_command(the_Command) @@ -137,7 +137,7 @@ def DEF(config): # must run with rank 1 processes = konfig["NUMBER_PART"] - the_Command = "SU2_DEF%s %s" % (quote, tempname) + the_Command = "SU2_DEF %s" % tempname the_Command = build_command(the_Command, processes) run_command(the_Command) @@ -164,7 +164,7 @@ def DOT(config): processes = konfig["NUMBER_PART"] - the_Command = "SU2_DOT_AD%s %s" % (quote, tempname) + the_Command = "SU2_DOT_AD %s" % tempname else: tempname = "config_DOT.cfg" @@ -172,7 +172,7 @@ def DOT(config): processes = konfig["NUMBER_PART"] - the_Command = "SU2_DOT%s %s" % (quote, tempname) + the_Command = "SU2_DOT %s" % tempname the_Command = build_command(the_Command, processes) run_command(the_Command) @@ -195,7 +195,7 @@ def GEO(config): # must run with rank 1 processes = konfig["NUMBER_PART"] - the_Command = "SU2_GEO%s %s" % (quote, tempname) + the_Command = "SU2_GEO %s" % tempname the_Command = build_command(the_Command, processes) run_command(the_Command) @@ -217,7 +217,7 @@ def SOL(config): # must run with rank 1 processes = konfig["NUMBER_PART"] - the_Command = "SU2_SOL%s %s" % (quote, tempname) + the_Command = "SU2_SOL %s" % tempname the_Command = build_command(the_Command, processes) run_command(the_Command) @@ -239,7 +239,7 @@ def SOL_FSI(config): # must run with rank 1 processes = konfig["NUMBER_PART"] - the_Command = "SU2_SOL%s %s 2" % (quote, tempname) + the_Command = "SU2_SOL %s 2" % tempname the_Command = build_command(the_Command, processes) run_command(the_Command) @@ -255,7 +255,14 @@ def SOL_FSI(config): def build_command(the_Command, processes=0): """builds an mpi command for given number of processes""" - the_Command = quote + (base_Command % the_Command) + # quote the executable path, SU2_RUN may contain spaces + executable, _, arguments = the_Command.partition(" ") + executable = base_Command % executable + if sys.platform == "win32": + executable = '"%s"' % executable + else: + executable = shlex.quote(executable) + the_Command = "%s %s" % (executable, arguments) if arguments else executable if processes > 1: if not mpi_Command: raise RuntimeError("could not find an mpi interface") diff --git a/config_template.cfg b/config_template.cfg index 5e611f22ed48..cb6647cd03d7 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1035,6 +1035,12 @@ FLAME_INIT_METHOD= FLAME_FRONT % respectively are included in the manifold. % By default, this option is disabled. PREFERENTIAL_DIFFUSION= NO +% +% Major species carrying the resolved Eq. (14) preferential diffusion fluxes of the +% SOURCE_TERM method. The manifold variable names are composed from these, so order and +% spelling must match the table: "H2" implies "D__H2" for every controlling variable +% , and the mass fraction "Y-H2". Required with PREFERENTIAL_DIFFUSION_METHOD= SOURCE_TERM. +PREFERENTIAL_DIFFUSION_MAJOR_SPECIES= (H2, H2O, H) % FLAME_FRONT initialization % the flame is initialized using a plane, defined by a point and a normal. On one side, the solution is initialized