diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 73b6b81e29e4..f34726036935 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1525,14 +1525,23 @@ static const MapType Species_Model_Map = { }; /*! - * \brief Progress variable and enthalpy are the first and second entries in the lookup table. - * \note The order matters. + * \brief Names of the flamelet controlling variables in CONTROLLING_VARIABLE_NAMES; the name sets the role. + */ +namespace FLAMELET_CV_NAME { + constexpr const char* PROGVAR = "ProgressVariable"; + constexpr const char* ENTHALPY = "EnthalpyTot"; + constexpr const char* MIXFRAC = "MixtureFraction"; + constexpr const char* PROGVAR_VAR = "ProgressVariableVariance"; + constexpr const char* MIXFRAC_VAR = "MixtureFractionVariance"; +} + +/*! + * \brief Manifold variables giving the range of the progress variable, the support of its PDF. */ -enum FLAMELET_SCALAR_VARIABLES { - I_PROGVAR, - I_ENTH, - I_MIXFRAC, -}; +namespace FLAMELET_PROGVAR_BOUNDS { + constexpr const char* MIN = "ProgressVariableMin"; + constexpr const char* MAX = "ProgressVariableMax"; +} /*! * \brief the source terms for the flamelet method. At the moment only progress variable @@ -1549,6 +1558,7 @@ enum FLAMELET_LOOKUP_OPS { PREFDIF, /*!< \brief Preferential diffusion scalars. */ SOURCES, /*!< \brief Scalar source terms (controlling variables, passive species).*/ LOOKUP, /*!< \brief Passive look-up variables specified in config. */ + PROGVAR_BOUNDS, /*!< \brief Minimum and maximum of the progress variable. */ }; /*! @@ -1607,6 +1617,12 @@ struct FluidFlamelet_ParsedOptions { unsigned short n_user_scalars = 0; /*!< \brief Number of user defined (auxiliary) scalar transport equations. */ unsigned short n_user_sources = 0; /*!< \brief Number of source terms for user defined (auxiliary) scalar transport equations. */ unsigned short n_control_vars = 0; /*!< \brief Number of controlling variables (independent variables) for the LUT. */ + int idx_progvar = -1; /*!< \brief Scalar index of the progress variable, -1 if absent. */ + int idx_enth = -1; /*!< \brief Scalar index of the total enthalpy. */ + int idx_mixfrac = -1; /*!< \brief Scalar index of the mixture fraction, -1 if absent. */ + int idx_progvar_var = -1; /*!< \brief Scalar index of the progress variable variance, -1 if absent. */ + int idx_mixfrac_var = -1; /*!< \brief Scalar index of the mixture fraction variance, -1 if absent. */ + su2double variance_dissipation_coeff = 2.0; /*!< \brief C_d in the variance dissipation C_d rho (eps/k) var. */ std::string *controlling_variable_names; /*!< \brief Names of the independent, transported scalars. */ std::string* cv_source_names; /*!< \brief Names of the source terms of the independent, transported scalars. */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index f2b171c9571f..635614c30bb3 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1528,6 +1528,9 @@ void CConfig::SetConfig_Options() { /*!\brief FLAME_INIT_IGNITION \n DESCRIPTION: Ignition temperature for the flame initialization \ingroup Config*/ addDoubleOption("FLAME_INIT_IGNITION", flamelet_ParsedOptions.Flame_T_ignition, 5000.0); + /*!\brief FLAMELET_VARIANCE_DISSIPATION_COEFF \n DESCRIPTION: Coefficient C_d of the scalar variance dissipation C_d rho (eps/k) var \ingroup Config*/ + addDoubleOption("FLAMELET_VARIANCE_DISSIPATION_COEFF", flamelet_ParsedOptions.variance_dissipation_coeff, 2.0); + /*--- Options related to mass diffusivity and thereby the species solver. ---*/ /*!\brief DIFFUSIVITY_MODEL\n DESCRIPTION: mass diffusivity model \n DEFAULT constant disffusivity \ingroup Config*/ @@ -6109,6 +6112,51 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i /*--- We can have additional user defined transported scalars ---*/ flamelet_ParsedOptions.n_scalars = flamelet_ParsedOptions.n_control_vars + flamelet_ParsedOptions.n_user_scalars; + /*--- The name of each controlling variable sets its role. ---*/ + auto& opts = flamelet_ParsedOptions; + for (int iCV = 0; iCV < opts.n_control_vars; iCV++) { + const string& name = opts.controlling_variable_names[iCV]; + int* idx = nullptr; + if (name == FLAMELET_CV_NAME::PROGVAR) idx = &opts.idx_progvar; + else if (name == FLAMELET_CV_NAME::ENTHALPY) idx = &opts.idx_enth; + else if (name == FLAMELET_CV_NAME::MIXFRAC) idx = &opts.idx_mixfrac; + else if (name == FLAMELET_CV_NAME::PROGVAR_VAR) idx = &opts.idx_progvar_var; + else if (name == FLAMELET_CV_NAME::MIXFRAC_VAR) idx = &opts.idx_mixfrac_var; + else + SU2_MPI::Error("Unknown controlling variable " + name + ", expected " + FLAMELET_CV_NAME::PROGVAR + ", " + + FLAMELET_CV_NAME::ENTHALPY + ", " + FLAMELET_CV_NAME::MIXFRAC + ", " + + FLAMELET_CV_NAME::PROGVAR_VAR + " or " + FLAMELET_CV_NAME::MIXFRAC_VAR + ".", + CURRENT_FUNCTION); + if (*idx >= 0) + SU2_MPI::Error("Controlling variable " + name + " is listed twice.", + CURRENT_FUNCTION); + *idx = iCV; + } + if (opts.idx_enth < 0) + SU2_MPI::Error(string("CONTROLLING_VARIABLE_NAMES must include ") + FLAMELET_CV_NAME::ENTHALPY + ".", + CURRENT_FUNCTION); + if (opts.idx_progvar < 0 && opts.idx_mixfrac < 0) + SU2_MPI::Error(string("CONTROLLING_VARIABLE_NAMES must include ") + FLAMELET_CV_NAME::PROGVAR + " or " + + FLAMELET_CV_NAME::MIXFRAC + ".", CURRENT_FUNCTION); + + /*--- The thickened flame correction measures the progress variable gradient. ---*/ + if (opts.idx_progvar < 0) opts.thickenedflame_correction = false; + + /*--- A variance needs its mean and the turbulence time scale omega of the SST model. ---*/ + if ((opts.idx_progvar_var >= 0) && (opts.idx_progvar < 0)) + SU2_MPI::Error(string(FLAMELET_CV_NAME::PROGVAR_VAR) + " requires " + FLAMELET_CV_NAME::PROGVAR + ".", + CURRENT_FUNCTION); + if ((opts.idx_mixfrac_var >= 0) && (opts.idx_mixfrac < 0)) + SU2_MPI::Error(string(FLAMELET_CV_NAME::MIXFRAC_VAR) + " requires " + FLAMELET_CV_NAME::MIXFRAC + ".", + CURRENT_FUNCTION); + if ((opts.idx_progvar_var >= 0) || (opts.idx_mixfrac_var >= 0)) { + if (Kind_Turb_Model != TURB_MODEL::SST) + SU2_MPI::Error("Scalar variance transport requires KIND_TURB_MODEL= SST.", CURRENT_FUNCTION); + if (opts.preferential_diffusion) + SU2_MPI::Error("PREFERENTIAL_DIFFUSION is not available with scalar variance transport.", CURRENT_FUNCTION); + opts.thickenedflame_correction = false; + } + /*--- Check that spark ignition has required parameters defined ---*/ if (flamelet_ParsedOptions.ignition_method == FLAMELET_INIT_TYPE::SPARK) { /*--- Check if SPARK_INIT was explicitly set in config file ---*/ diff --git a/SU2_CFD/include/fluid/CFluidFlamelet.hpp b/SU2_CFD/include/fluid/CFluidFlamelet.hpp index 88b2e6603ee6..664edfa6ffcc 100644 --- a/SU2_CFD/include/fluid/CFluidFlamelet.hpp +++ b/SU2_CFD/include/fluid/CFluidFlamelet.hpp @@ -48,8 +48,7 @@ class CFluidFlamelet final : public CFluidModel { su2double beta_progvar, beta_enth_thermal, beta_enth, beta_mixfrac; int rank; - bool include_mixture_fraction = false, /*!< \brief include mixture fraction in controlling variables. */ - preferential_diffusion = false; /*!< \brief use preferential diffusion physics. */ + bool preferential_diffusion = false; /*!< \brief use preferential diffusion physics. */ unsigned short n_scalars, n_lookups, n_user_scalars, /*!< \brief number of passive reactant species. */ n_control_vars; /*!< \brief number of controlling variables. */ @@ -68,7 +67,12 @@ class CFluidFlamelet final : public CFluidModel { vector LUT_idx_TD, LUT_idx_Sources, LUT_idx_LookUp, - LUT_idx_PD; + LUT_idx_PD, + LUT_idx_ProgVarBounds; + + vector manifold_inputs; /*!< \brief Controlling variables with the variances normalized. */ + su2double limit_progvar_var = 0, /*!< \brief Upper limit of the progress variable variance. */ + limit_mixfrac_var = 0; /*!< \brief Upper limit of the mixture fraction variance. */ /*--- Class variables for the multi-layer perceptron method ---*/ #ifdef USE_MLPCPP @@ -78,7 +82,7 @@ class CFluidFlamelet final : public CFluidModel { MLPToolbox::CIOMap* iomap_PD; /*!< \brief Input-output map for the preferential diffusion scalars. */ MLPToolbox::CIOMap* iomap_Sources; /*!< \brief Input-output map for species source terms. */ MLPToolbox::CIOMap* iomap_LookUp; /*!< \brief Input-output map for passive look-up terms. */ - MLPToolbox::CIOMap* iomap_Current; + MLPToolbox::CIOMap* iomap_ProgVarBounds; /*!< \brief Input-output map for the progress variable range. */ #endif vector scalars_vector; @@ -86,15 +90,33 @@ class CFluidFlamelet final : public CFluidModel { vector varnames_TD, /*!< \brief Lookup names for thermodynamic state variables. */ varnames_Sources, /*!< \brief Lookup names for source terms. */ varnames_LookUp, /*!< \brief Lookup names for passive look-up terms. */ - varnames_PD; /*!< \brief Lookup names for preferential diffusion scalars. */ + varnames_PD, /*!< \brief Lookup names for preferential diffusion scalars. */ + varnames_ProgVarBounds; /*!< \brief Lookup names for the minimum and maximum progress variable. */ vector val_vars_TD, /*!< \brief References to thermodynamic state variables. */ val_vars_Sources, /*!< \brief References to source terms. */ val_vars_LookUp, /*!< \brief References passive look-up terms. */ - val_vars_PD; /*!< \brief References to preferential diffusion scalars. */ + val_vars_PD, /*!< \brief References to preferential diffusion scalars. */ + val_vars_ProgVarBounds; /*!< \brief References to the minimum and maximum progress variable. */ void PreprocessLookUp(CConfig* config); + /*! + * \brief Interpolate the manifold, without normalizing the inputs. + * \param[in] lookup_type - Look-up operation to be performed (FLAMELET_LOOKUP_OPS). + * \param[in] inputs - Manifold inputs in the order of the controlling variables. + * \param[out] output_refs - Interpolated output variables. + * \returns Whether the query is inside the manifold. + */ + bool LookUpManifold(unsigned short lookup_type, const vector& inputs, vector& output_refs); + + /*! + * \brief Set manifold_inputs to the controlling variables, each variance divided by its upper limit. + * \param[in] scalars - Transported scalars. + * \returns Whether the look-up of the progress variable range is inside the manifold. + */ + bool SetManifoldInputs(const su2double* scalars); + /*! \brief * Returns true if the string is null or zero (ignores case). */ @@ -135,6 +157,14 @@ class CFluidFlamelet final : public CFluidModel { unsigned long EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, vector& output_refs); + /*! + * \brief Largest value a scalar variance can take, (mean - min)(max - mean) for its mean. + * \param[in] scalars - Transported scalars. + * \param[in] i_var - Index of the variance. + * \returns Upper limit of the variance. + */ + su2double GetScalarVarianceLimit(const su2double* scalars, unsigned short i_var) override; + /*! * \brief Check for out-of-bounds condition for data set interpolation. * \return - within bounds (0) or out of bounds (1). diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index af9aa3e059b7..47528e25ff39 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/SU2_CFD/include/fluid/CFluidModel.hpp @@ -161,6 +161,14 @@ class CFluidModel { return 0; } + /*! + * \brief Largest value a scalar variance can take, (mean - min)(max - mean) for its mean. + * \param[in] scalars - Transported scalars. + * \param[in] i_var - Index of the variance. + * \returns Upper limit of the variance. + */ + virtual su2double GetScalarVarianceLimit(const su2double* scalars, unsigned short i_var) { return 1.0e15; } + /*! * \brief Get fluid dynamic viscosity. */ diff --git a/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp b/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp index dba70c948816..0d58cf2df681 100644 --- a/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp +++ b/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp @@ -50,7 +50,9 @@ 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), + idxProgVar(config.GetFlameletParsedOptions().idx_progvar), + idxEnth(config.GetFlameletParsedOptions().idx_enth) {} /*! * \brief Preferential diffusion, two terms with the shape of the ordinary diffusion but of @@ -114,8 +116,6 @@ class CScalarFlux_Flamelet final } /*--- Thermal term, on the enthalpy equation alone, driven by the temperature gradient. ---*/ - if (nControlVars <= I_ENTH) return; - const Double T_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.Temperature()); const Double T_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.Temperature()); @@ -123,33 +123,29 @@ class CScalarFlux_Flamelet final const auto gradT_j = gatherVariables(jPoint, side_j.flowNodes->GetGradient_Primitive(), idx.Temperature()); const Double Dth_i = gatherVariables(iPoint, side_i.scalarNodes.GetAuxVar(), I_BETA_ENTH_THERMAL) * - gatherVariables(iPoint, side_i.scalarNodes.GetDiffusivity(), I_ENTH); + gatherVariables(iPoint, side_i.scalarNodes.GetDiffusivity(), idxEnth); const Double Dth_j = gatherVariables(jPoint, side_j.scalarNodes.GetAuxVar(), I_BETA_ENTH_THERMAL) * - gatherVariables(jPoint, side_j.scalarNodes.GetDiffusivity(), I_ENTH); + gatherVariables(jPoint, side_j.scalarNodes.GetDiffusivity(), idxEnth); const Double Dth = 0.5 * (rho.i * Dth_i + rho.j * Dth_j) + diffTurb; const Double projGradT = projectedGradient(opt, gradT_i, gradT_j, T_i, T_j, normal, vector_ij, dist2_ij); - res.flux_i(I_ENTH) -= Dth * projGradT; - if (!opt.oneSided) res.flux_j(I_ENTH) += Dth * projGradT; + res.flux_i(idxEnth) -= Dth * projGradT; + if (!opt.oneSided) res.flux_j(idxEnth) += Dth * projGradT; } private: const bool preferentialDiffusion; const unsigned short nControlVars; + const int idxProgVar, idxEnth; /*! * \brief Auxiliary variable holding the beta scalar of a controlling variable. */ - static FORCEINLINE unsigned short betaIndex(unsigned short iScalar) { - switch (iScalar) { - case I_PROGVAR: - return I_BETA_PROGVAR; - case I_ENTH: - return I_BETA_ENTH; - default: - return I_BETA_MIXFRAC; - } + FORCEINLINE unsigned short betaIndex(int iScalar) const { + if (iScalar == idxProgVar) return I_BETA_PROGVAR; + if (iScalar == idxEnth) return I_BETA_ENTH; + return I_BETA_MIXFRAC; } /*! diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index 7543d6d7f282..80a75b6b26f2 100644 --- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp @@ -41,7 +41,10 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { const su2double default_flame_thickness{1.0}; su2double global_flame_thickness; bool calc_flame_thickness{false}; - bool include_mixture_fraction = false; /*!< \brief include mixture fraction as a controlling variable. */ + int idx_progvar = -1; /*!< \brief Scalar index of the progress variable, -1 if absent. */ + int idx_enth = -1; /*!< \brief Scalar index of the total enthalpy. */ + int idx_flame_front = -1; /*!< \brief Scalar set by FLAME_FRONT: progress variable, else mixture fraction. */ + vector> variance_sources; /*!< \brief Variance and mean indices of variances with a chemical source. */ /*!< \brief Number of points outside the manifold domain, shared across OpenMP threads so it can be * accumulated atomically and reduced by the master thread alone in Preprocessing. */ unsigned long n_not_in_domain_local = 0; @@ -51,8 +54,27 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ + /*! + * \brief Limit the scalar variances to the maximum of a two-point distribution on the range of their mean. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Postprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, + unsigned short iMesh) override; + void SetPreconditioner(CGeometry* geometry, CSolver** solver_container, CConfig* config); + /*! + * \brief Add the turbulent production 2 mu_t/Sc_t |grad(mean)|^2 and dissipation C_d rho beta* omega var + * of a scalar variance to the residual, the dissipation implicitly. + * \param[in] i_var - Scalar index of the variance. + * \param[in] i_mean - Scalar index of its mean. + */ + void ScalarVarianceResidual(CGeometry* geometry, CSolver** solver_container, const CConfig* config, int i_var, + int i_mean); + /*! * \brief Generic implementation of the isothermal wall also covering CHT cases, * for which the wall temperature is given by GetConjugateHeatVariable. diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index f142da539175..1e959fcae847 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -44,7 +44,6 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati /* -- number of auxiliary species transport equations, e.g. 1=CO, 2=NOx --- */ n_user_scalars = flamelet_options.n_user_scalars; n_control_vars = flamelet_options.n_control_vars; - include_mixture_fraction = (n_control_vars == 3); n_scalars = flamelet_options.n_scalars; if (rank == MASTER_NODE) { @@ -53,6 +52,7 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati cout << "Number of control variables: " << n_control_vars << endl; } scalars_vector.resize(n_scalars); + manifold_inputs.resize(n_control_vars); table_scalar_names.resize(n_scalars); for (auto iCV = 0u; iCV < n_control_vars; iCV++) @@ -73,13 +73,15 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: + if (n_control_vars > 3) + SU2_MPI::Error("The look-up table supports at most three controlling variables.", CURRENT_FUNCTION); if (rank == MASTER_NODE) { cout << "*****************************************" << endl; cout << "*** initializing the lookup table ***" << endl; cout << "*****************************************" << endl; } - look_up_table = new CLookUpTable(datadriven_fluid_options.datadriven_filenames[0], table_scalar_names[I_PROGVAR], - table_scalar_names[I_ENTH]); + look_up_table = new CLookUpTable(datadriven_fluid_options.datadriven_filenames[0], table_scalar_names[0], + table_scalar_names[1]); break; default: @@ -116,13 +118,14 @@ CFluidFlamelet::~CFluidFlamelet() { delete iomap_LookUp; delete lookup_mlp; if (preferential_diffusion) delete iomap_PD; + if (flamelet_options.idx_progvar_var >= 0) delete iomap_ProgVarBounds; } #endif } void CFluidFlamelet::SetTDState_h(su2double val_enthalpy, const su2double* val_scalars) { /*--- For the fluid flamelet model, the enthalpy (and the temperature) are passive scalars. - val_scalars contains the enthalpy as the second variable: val_scalars= (Progress_variable, enthalpy,...). + val_scalars contains the enthalpy at the index of EnthalpyTot in CONTROLLING_VARIABLE_NAMES. Consequently, the energy equation is not solved when the fluid flamelet model is used; instead the energy equation for enthalpy is solved in the species flamelet solver, and the enthalpy solution is overwritten in the CIncEulerSolver. Likewise, The temperature is retrieved from the look up table. Then, the thermodynamics state is fully determined with @@ -137,7 +140,7 @@ void CFluidFlamelet::SetTDState_T(su2double val_temperature, const su2double* va /*--- Add all quantities and their names to the look up vectors. ---*/ EvaluateDataSet(scalars_vector, FLAMELET_LOOKUP_OPS::THERMO, val_vars_TD); - Enthalpy = scalars_vector[1]; + Enthalpy = scalars_vector[flamelet_options.idx_enth]; Temperature = val_vars_TD[LOOKUP_TD::TEMPERATURE]; Cp = val_vars_TD[LOOKUP_TD::HEATCAPACITY]; Mu = val_vars_TD[LOOKUP_TD::VISCOSITY]; @@ -216,16 +219,23 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { varnames_PD.resize(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); val_vars_PD.resize(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); - varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR] = "Beta_ProgVar"; + /*--- Beta scalars of absent controlling variables are not in the manifold. ---*/ + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR] = (flamelet_options.idx_progvar >= 0) ? "Beta_ProgVar" : "NULL"; 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_MIXFRAC] = (flamelet_options.idx_mixfrac >= 0) ? "Beta_MixFrac" : "NULL"; 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; + /*--- The progress variable range normalizes its variance. ---*/ + if (flamelet_options.idx_progvar_var >= 0) { + varnames_ProgVarBounds = {FLAMELET_PROGVAR_BOUNDS::MIN, FLAMELET_PROGVAR_BOUNDS::MAX}; + val_vars_ProgVarBounds.resize(varnames_ProgVarBounds.size()); + } + preferential_diffusion = flamelet_options.preferential_diffusion; if (!preferential_diffusion && flamelet_options.preferential_diffusion) @@ -238,12 +248,16 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { iomap_LookUp = new MLPToolbox::CIOMap(controlling_variable_names, varnames_LookUp); lookup_mlp->PairVariableswithMLPs(*iomap_TD); lookup_mlp->PairVariableswithMLPs(*iomap_Sources); - if (n_lookups > 1) + if (n_lookups > 0) lookup_mlp->PairVariableswithMLPs(*iomap_LookUp); if (preferential_diffusion) { iomap_PD = new MLPToolbox::CIOMap(controlling_variable_names, varnames_PD); lookup_mlp->PairVariableswithMLPs(*iomap_PD); } + if (!varnames_ProgVarBounds.empty()) { + iomap_ProgVarBounds = new MLPToolbox::CIOMap(controlling_variable_names, varnames_ProgVarBounds); + lookup_mlp->PairVariableswithMLPs(*iomap_ProgVarBounds); + } #endif } else { for (auto iVar = 0u; iVar < varnames_TD.size(); iVar++) { @@ -268,76 +282,131 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { } if (preferential_diffusion) { for (auto iVar = 0u; iVar < varnames_PD.size(); iVar++) { - LUT_idx_PD.push_back(look_up_table->GetIndexOfVar(varnames_PD[iVar])); + LUT_idx_PD.push_back(noSource(varnames_PD[iVar]) ? look_up_table->GetNullIndex() + : look_up_table->GetIndexOfVar(varnames_PD[iVar])); } } + if (!varnames_ProgVarBounds.empty()) { + if (!look_up_table->CheckForVariables(varnames_ProgVarBounds)) + SU2_MPI::Error(string("Transporting ") + FLAMELET_CV_NAME::PROGVAR_VAR + " requires the variables " + + FLAMELET_PROGVAR_BOUNDS::MIN + " and " + FLAMELET_PROGVAR_BOUNDS::MAX + " in the table.", + CURRENT_FUNCTION); + for (const auto& name : varnames_ProgVarBounds) LUT_idx_ProgVarBounds.push_back(look_up_table->GetIndexOfVar(name)); + } } } -unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, - vector& output_refs) { - AD::StartPreacc(); - for (auto iVar = 0u; iVar < input_scalar.size(); iVar++) AD::SetPreaccIn(input_scalar[iVar]); +su2double CFluidFlamelet::GetScalarVarianceLimit(const su2double* scalars, unsigned short i_var) { + SetManifoldInputs(scalars); + if (i_var == flamelet_options.idx_progvar_var) return limit_progvar_var; + if (i_var == flamelet_options.idx_mixfrac_var) return limit_mixfrac_var; + return CFluidModel::GetScalarVarianceLimit(scalars, i_var); +} - su2double val_enth = input_scalar[I_ENTH]; - su2double val_prog = input_scalar[I_PROGVAR]; - su2double val_mixfrac = include_mixture_fraction ? input_scalar[I_MIXFRAC] : 0.0; - vector val_vars; - vector refs_vars; - vector LUT_idx; +bool CFluidFlamelet::SetManifoldInputs(const su2double* scalars) { + for (auto iCV = 0u; iCV < n_control_vars; iCV++) manifold_inputs[iCV] = scalars[iCV]; + + /*--- Variance divided by its upper limit, 0 where the mean is at a bound. ---*/ + auto normalize = [](su2double variance, su2double limit) { + return (limit > 0) ? fmin(fmax(variance / limit, su2double(0.0)), su2double(1.0)) : su2double(0.0); + }; + + const int i_zvar = flamelet_options.idx_mixfrac_var; + if (i_zvar >= 0) { + const su2double z = scalars[flamelet_options.idx_mixfrac]; + limit_mixfrac_var = fmax(z * (1.0 - z), su2double(0.0)); + manifold_inputs[i_zvar] = normalize(scalars[i_zvar], limit_mixfrac_var); + } + + bool inside{true}; + const int i_cvar = flamelet_options.idx_progvar_var; + if (i_cvar >= 0) { + /*--- The range of the progress variable does not depend on its variance. ---*/ + manifold_inputs[i_cvar] = 0.0; + inside = LookUpManifold(FLAMELET_LOOKUP_OPS::PROGVAR_BOUNDS, manifold_inputs, val_vars_ProgVarBounds); + const su2double c = scalars[flamelet_options.idx_progvar]; + limit_progvar_var = fmax((c - val_vars_ProgVarBounds[0]) * (val_vars_ProgVarBounds[1] - c), su2double(0.0)); + manifold_inputs[i_cvar] = normalize(scalars[i_cvar], limit_progvar_var); + } + return inside; +} + +bool CFluidFlamelet::LookUpManifold(unsigned short lookup_type, const vector& inputs, + vector& output_refs) { + vector* LUT_idx = nullptr; +#ifdef USE_MLPCPP + MLPToolbox::CIOMap* iomap = nullptr; +#endif switch (lookup_type) { case FLAMELET_LOOKUP_OPS::THERMO: - LUT_idx = LUT_idx_TD; + LUT_idx = &LUT_idx_TD; #ifdef USE_MLPCPP - iomap_Current = iomap_TD; + iomap = iomap_TD; #endif break; case FLAMELET_LOOKUP_OPS::PREFDIF: - LUT_idx = LUT_idx_PD; + LUT_idx = &LUT_idx_PD; #ifdef USE_MLPCPP - iomap_Current = iomap_PD; + iomap = iomap_PD; #endif break; case FLAMELET_LOOKUP_OPS::SOURCES: - LUT_idx = LUT_idx_Sources; + LUT_idx = &LUT_idx_Sources; #ifdef USE_MLPCPP - iomap_Current = iomap_Sources; + iomap = iomap_Sources; #endif break; case FLAMELET_LOOKUP_OPS::LOOKUP: - LUT_idx = LUT_idx_LookUp; + LUT_idx = &LUT_idx_LookUp; #ifdef USE_MLPCPP - iomap_Current = iomap_LookUp; + iomap = iomap_LookUp; #endif break; - default: + case FLAMELET_LOOKUP_OPS::PROGVAR_BOUNDS: + LUT_idx = &LUT_idx_ProgVarBounds; +#ifdef USE_MLPCPP + iomap = iomap_ProgVarBounds; +#endif break; + default: + SU2_MPI::Error("Unknown manifold look-up operation.", CURRENT_FUNCTION); } - /*--- Add all quantities and their names to the look up vectors. ---*/ bool inside{true}; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: - if (output_refs.size() != LUT_idx.size()) + if (output_refs.size() != LUT_idx->size()) SU2_MPI::Error(string("Output vector size incompatible with manifold lookup operation."), CURRENT_FUNCTION); - if (include_mixture_fraction) { - inside = look_up_table->LookUp_XYZ(LUT_idx, output_refs, val_prog, val_enth, val_mixfrac); + /*--- The table axes follow the order of the controlling variables, the third one is the table level. ---*/ + if (n_control_vars == 3) { + inside = look_up_table->LookUp_XYZ(*LUT_idx, output_refs, inputs[0], inputs[1], inputs[2]); } else { - inside = look_up_table->LookUp_XY(LUT_idx, output_refs, val_prog, val_enth); + inside = look_up_table->LookUp_XY(*LUT_idx, output_refs, inputs[0], inputs[1]); } break; - case ENUM_DATADRIVEN_METHOD::MLP: - refs_vars.resize(output_refs.size()); - for (auto iVar = 0u; iVar < output_refs.size(); iVar++) refs_vars[iVar] = &output_refs[iVar]; + case ENUM_DATADRIVEN_METHOD::MLP: { #ifdef USE_MLPCPP - inside=lookup_mlp->Predict(*iomap_Current, input_scalar, refs_vars); + vector refs_vars(output_refs.size()); + for (auto iVar = 0u; iVar < output_refs.size(); iVar++) refs_vars[iVar] = &output_refs[iVar]; + inside = lookup_mlp->Predict(*iomap, inputs, refs_vars); #endif break; + } default: break; } - if (inside) extrapolation = 0; - else extrapolation = 1; + return inside; +} + +unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, + vector& output_refs) { + AD::StartPreacc(); + for (auto iVar = 0u; iVar < input_scalar.size(); iVar++) AD::SetPreaccIn(input_scalar[iVar]); + + bool inside = SetManifoldInputs(input_scalar.data()); + inside = LookUpManifold(lookup_type, manifold_inputs, output_refs) && inside; + + extrapolation = inside ? 0 : 1; for (auto iVar = 0u; iVar < output_refs.size(); iVar++) AD::SetPreaccOut(output_refs[iVar]); AD::EndPreacc(); return extrapolation; diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 19f07df2e598..b36fac00d48f 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -3948,9 +3948,12 @@ void CSolver::LoadInletProfile(CGeometry **geometry, break; case SPECIES_MODEL::FLAMELET: { const auto& flamelet_config_options = config->GetFlameletParsedOptions(); - /*--- 2-equation flamelet model ---*/ - columnName << left << setw(24) << "PROGRESSVAR" << left << setw(24) << "ENTHALPYTOT"; - columnValue << config->GetInlet_SpeciesVal(Marker_Tag)[0] << "\t" << config->GetInlet_SpeciesVal(Marker_Tag)[1] <<"\t"; + for (unsigned short iCV = 0; iCV < flamelet_config_options.n_control_vars; iCV++) { + string cv_name = flamelet_config_options.controlling_variable_names[iCV]; + transform(cv_name.begin(), cv_name.end(), cv_name.begin(), ::toupper); + columnName << left << setw(24) << cv_name; + columnValue << config->GetInlet_SpeciesVal(Marker_Tag)[iCV] << "\t"; + } /*--- auxiliary species transport equations ---*/ for (unsigned short iReactant = 0; iReactant < flamelet_config_options.n_user_scalars; iReactant++) { columnName << left << setw(24) << flamelet_config_options.user_scalar_names[iReactant]; diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 9fdf81663b4e..d3499c0c723b 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -48,10 +48,24 @@ CSpeciesFlameletSolver::CSpeciesFlameletSolver(CGeometry* geometry, CConfig* con /*--- Dimension of the problem. ---*/ nVar = flamelet_config_options.n_scalars; - include_mixture_fraction = (flamelet_config_options.n_control_vars == 3); + idx_progvar = flamelet_config_options.idx_progvar; + idx_enth = flamelet_config_options.idx_enth; + idx_flame_front = (idx_progvar >= 0) ? idx_progvar : flamelet_config_options.idx_mixfrac; Initialize(geometry, config, iMesh, nVar); + /*--- A variance cannot become negative. ---*/ + for (const int i_var : {flamelet_config_options.idx_progvar_var, flamelet_config_options.idx_mixfrac_var}) + if (i_var >= 0) lowerlimit[i_var] = max(lowerlimit[i_var], su2double(0.0)); + + for (const auto& var_mean : {make_pair(flamelet_config_options.idx_progvar_var, idx_progvar), + make_pair(flamelet_config_options.idx_mixfrac_var, flamelet_config_options.idx_mixfrac)}) { + if (var_mean.first < 0) continue; + string source_name = flamelet_config_options.cv_source_names[var_mean.first]; + transform(source_name.begin(), source_name.end(), source_name.begin(), ::toupper); + if (source_name != "NULL" && source_name != "ZERO") variance_sources.push_back(var_mean); + } + /*--- Initialize the solution to the far-field state everywhere. ---*/ nodes = new CSpeciesFlameletVariable(Solution_Inf, nPoint, nDim, nVar, config); @@ -214,7 +228,7 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** } const su2double temp_inlet = config->GetInc_Temperature_Init(); - su2double enth_inlet = config->GetSpecies_Init()[I_ENTH]; + su2double enth_inlet = config->GetSpecies_Init()[idx_enth]; su2double prog_burnt = 0, prog_unburnt, point_loc; su2double scalar_init[MAXNVAR]= {0.0}; @@ -254,11 +268,11 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** /*--- Set enthalpy based on initial temperature and scalars. ---*/ n_not_iterated_local += GetEnthFromTemp(fluid_model_local, temp_inlet, config->GetSpecies_Init(), &enth_inlet); - scalar_init[I_ENTH] = enth_inlet; + scalar_init[idx_enth] = enth_inlet; if (flame_front_ignition) prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init, flamelet_config_options.Flame_T_ignition); - prog_unburnt = config->GetSpecies_Init()[I_PROGVAR]; + prog_unburnt = config->GetSpecies_Init()[idx_flame_front]; const auto nPoint_iMesh = geometry[i_mesh]->GetnPoint(); SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long i_point = 0; i_point < nPoint_iMesh; i_point++) { @@ -277,26 +291,26 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** /* --- Unburnt region upstream of the flame. --- */ if (point_loc <= 0) { - scalar_init[I_PROGVAR] = prog_unburnt; + scalar_init[idx_flame_front] = prog_unburnt; n_points_unburnt_local++; /* --- Flame zone where we lineary increase from unburnt to burnt conditions. --- */ } else if ((point_loc > 0) && (point_loc <= flame_thickness)) { - scalar_init[I_PROGVAR] = prog_unburnt + point_loc * (prog_burnt - prog_unburnt) / flame_thickness; + scalar_init[idx_flame_front] = prog_unburnt + point_loc * (prog_burnt - prog_unburnt) / flame_thickness; n_points_flame_local++; /* --- Burnt region behind the flame zone. --- */ } else if ((point_loc > flame_thickness) && (point_loc <= flame_thickness + flame_burnt_thickness)) { - scalar_init[I_PROGVAR] = prog_burnt; + scalar_init[idx_flame_front] = prog_burnt; n_points_burnt_local++; /* --- Unburnt region downstream of the flame and burnt region. --- */ } else { - scalar_init[I_PROGVAR] = prog_unburnt; + scalar_init[idx_flame_front] = prog_unburnt; n_points_unburnt_local++; } } else { - scalar_init[I_PROGVAR] = prog_unburnt; + scalar_init[idx_flame_front] = prog_unburnt; } /* --- Perform manifold evaluation to check whether initial scalar is out-of-bounds. --- */ fluid_model_local->SetTDState_T(temp_inlet, scalar_init); @@ -419,6 +433,25 @@ void CSpeciesFlameletSolver::SetPreconditioner(CGeometry* geometry, CSolver** so END_SU2_OMP_FOR } +void CSpeciesFlameletSolver::Postprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, + unsigned short iMesh) { + SU2_ZONE_SCOPED + + const int variances[] = {flamelet_config_options.idx_progvar_var, flamelet_config_options.idx_mixfrac_var}; + if (variances[0] < 0 && variances[1] < 0) return; + + 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(); + for (const int i_var : variances) { + if (i_var < 0) continue; + const su2double limit = fluid_model_local->GetScalarVarianceLimit(nodes->GetSolution(i_point), i_var); + if (nodes->GetSolution(i_point, i_var) > limit) nodes->SetSolution(i_point, i_var, limit); + } + } + END_SU2_OMP_FOR +} + void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config, unsigned short iMesh) { SU2_ZONE_SCOPED @@ -427,10 +460,12 @@ void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solv const auto n_CV = flamelet_config_options.n_control_vars; const auto n_aux = flamelet_config_options.n_user_scalars; const auto* fn = static_cast(nodes); + const auto* flowNodes = solver_container[FLOW_SOL]->GetNodes(); SU2_OMP_FOR_STAT(omp_chunk_size) for (auto i_point = 0u; i_point < nPointDomain; i_point++) { const su2double volume = geometry->nodes->GetVolume(i_point); + const su2double density = flowNodes->GetDensity(i_point); /*--- Add source terms from the lookup table directly to the residual. ---*/ for (auto i_var = 0; i_var < nVar; i_var++) { @@ -439,21 +474,56 @@ void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solv /*--- Implicit: analytic Jacobian for auxiliary species from the split source form. * S_aux_i = source_prod_i + source_cons_i * Y_aux_i - * dS_aux_i/dY_aux_i = source_cons_i - * J_ii += -source_cons_i * V ---*/ + * dS_aux_i/d(rho Y_aux_i) = source_cons_i / rho + * J_ii += -source_cons_i / rho * V ---*/ if (implicit) { for (auto i_aux = 0u; i_aux < n_aux; i_aux++) { - Jacobian.AddVal2Diag(i_point, n_CV + i_aux, -fn->GetAuxSourceCons(i_point, i_aux) * volume); + Jacobian.AddVal2Diag(i_point, n_CV + i_aux, -fn->GetAuxSourceCons(i_point, i_aux) / density * volume); } } } END_SU2_OMP_FOR + if (flamelet_config_options.idx_progvar_var >= 0) + ScalarVarianceResidual(geometry, solver_container, config, flamelet_config_options.idx_progvar_var, idx_progvar); + if (flamelet_config_options.idx_mixfrac_var >= 0) + ScalarVarianceResidual(geometry, solver_container, config, flamelet_config_options.idx_mixfrac_var, + flamelet_config_options.idx_mixfrac); + /*--- call the species solver for the shared sources (axisymmetric and custom python source term) ---*/ CSpeciesSolver::Source_Residual(geometry, solver_container, numerics_container, config, iMesh); } +void CSpeciesFlameletSolver::ScalarVarianceResidual(CGeometry* geometry, CSolver** solver_container, + const CConfig* config, int i_var, int i_mean) { + SU2_ZONE_SCOPED + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const su2double Sc_t = config->GetSchmidt_Number_Turbulent(); + const su2double C_d = flamelet_config_options.variance_dissipation_coeff; + const su2double beta_star = solver_container[TURB_SOL]->GetConstants()[6]; + const auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + const auto* turbNodes = solver_container[TURB_SOL]->GetNodes(); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto i_point = 0u; i_point < nPointDomain; i_point++) { + const su2double volume = geometry->nodes->GetVolume(i_point); + + su2double grad_mean2 = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) grad_mean2 += pow(nodes->GetGradient(i_point, i_mean, iDim), 2); + const su2double production = 2.0 * flowNodes->GetEddyViscosity(i_point) / Sc_t * grad_mean2; + + /*--- Decay rate C_d eps/k with eps/k = beta* omega; the unknown is rho*var. ---*/ + const su2double decay = C_d * beta_star * turbNodes->GetSolution(i_point, 1); + const su2double rho_var = flowNodes->GetDensity(i_point) * nodes->GetSolution(i_point, i_var); + + LinSysRes(i_point, i_var) -= (production - decay * rho_var) * volume; + if (implicit) Jacobian.AddVal2Diag(i_point, i_var, decay * volume); + } + END_SU2_OMP_FOR +} + void CSpeciesFlameletSolver::BC_HeatFlux_Wall(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) { @@ -490,7 +560,7 @@ void CSpeciesFlameletSolver::BC_HeatFlux_Wall(CGeometry* geometry, CSolver** sol /*--- Neumann condition: q_wall is the prescribed heat flux (W/m^2, positive into domain). This adds a source term dH/dn * lambda = q_wall to the enthalpy residual. ---*/ - LinSysRes(iPoint, I_ENTH) -= Wall_HeatFlux * Area; + LinSysRes(iPoint, idx_enth) -= Wall_HeatFlux * Area; } END_SU2_OMP_FOR } @@ -511,12 +581,12 @@ void CSpeciesFlameletSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_cont } else { /*--- Use the enthalpy value directly from MARKER_INLET_SPECIES (default). The user is responsible for providing a thermodynamically consistent value. ---*/ - enth_inlet = config->GetInlet_SpeciesVal(Marker_Tag)[I_ENTH]; + enth_inlet = config->GetInlet_SpeciesVal(Marker_Tag)[idx_enth]; } SU2_OMP_FOR_STAT(OMP_MIN_SIZE) for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - Inlet_SpeciesVars[val_marker][iVertex][I_ENTH] = enth_inlet; + Inlet_SpeciesVars[val_marker][iVertex][idx_enth] = enth_inlet; } END_SU2_OMP_FOR @@ -550,7 +620,7 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo if (geometry->nodes->GetDomain(iPoint)) { if (config->GetMarker_StrongBC(Marker_Tag)) { /*--- Initial guess for enthalpy value. ---*/ - enth_wall = nodes->GetSolution(iPoint, I_ENTH); + enth_wall = nodes->GetSolution(iPoint, idx_enth); /*--- Set enthalpy on the wall. ---*/ n_not_iterated += GetEnthFromTemp(fluid_model_local, temp_wall, nodes->GetSolution(iPoint), &enth_wall); @@ -559,15 +629,15 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo condition (Dirichlet) and remove any contribution to the residual at this node. ---*/ - nodes->SetSolution(iPoint, I_ENTH, enth_wall); - nodes->SetSolution_Old(iPoint, I_ENTH, enth_wall); + nodes->SetSolution(iPoint, idx_enth, enth_wall); + nodes->SetSolution_Old(iPoint, idx_enth, enth_wall); - LinSysRes(iPoint, I_ENTH) = 0.0; + LinSysRes(iPoint, idx_enth) = 0.0; - nodes->SetVal_ResTruncError_Zero(iPoint, I_ENTH); + nodes->SetVal_ResTruncError_Zero(iPoint, idx_enth); if (implicit) { - Jacobian.DeleteValsRowi(iPoint, I_ENTH); + Jacobian.DeleteValsRowi(iPoint, idx_enth); } } else { /*--- Weak BC formulation. ---*/ @@ -597,7 +667,7 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo /*--- Apply a weak boundary condition for the energy equation. Compute the residual due to the prescribed heat flux. ---*/ - LinSysRes(iPoint, I_ENTH) -= thermal_conductivity * dTdn * Area; + LinSysRes(iPoint, idx_enth) -= thermal_conductivity * dTdn * Area; } } } @@ -641,10 +711,14 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF vector table_sources(flamelet_config_options.n_control_vars + 2 * flamelet_config_options.n_user_scalars); unsigned long misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::SOURCES, table_sources); - table_sources[I_PROGVAR] = fmax(0, table_sources[I_PROGVAR]); - nodes->SetTableMisses(iPoint, misses); - /*--- The source term for progress variable is always positive, we clip from below to make sure. --- */ + if (idx_progvar >= 0) table_sources[idx_progvar] = fmax(0, table_sources[idx_progvar]); + + /*--- Variance source 2 (mean(phi S) - mean(phi) mean(S)); the manifold gives mean(phi S). ---*/ + for (const auto& var_mean : variance_sources) + table_sources[var_mean.first] = + 2.0 * (table_sources[var_mean.first] - scalars[var_mean.second] * table_sources[var_mean.second]); + nodes->SetTableMisses(iPoint, misses); vector source_scalar(flamelet_config_options.n_scalars); for (auto iCV = 0u; iCV < flamelet_config_options.n_control_vars; iCV++) source_scalar[iCV] = table_sources[iCV]; @@ -715,7 +789,7 @@ unsigned long CSpeciesFlameletSolver::GetEnthFromTemp(CFluidModel* fluid_model, SU2_ZONE_SCOPED /*--- convergence criterion for temperature in [K], high accuracy needed for restarts. ---*/ su2double delta_temp_final = 0.001; - su2double enth_iter = scalar_solution[I_ENTH]; + su2double enth_iter = scalar_solution[idx_enth]; su2double delta_enth; su2double delta_temp_iter = 1e10; unsigned long exit_code = 0; @@ -730,7 +804,7 @@ unsigned long CSpeciesFlameletSolver::GetEnthFromTemp(CFluidModel* fluid_model, while ((abs(delta_temp_iter) > delta_temp_final) && (counter++ < counter_limit)) { /*--- Add all quantities and their names to the look up vectors. ---*/ - val_scalars[I_ENTH] = enth_iter; + val_scalars[idx_enth] = enth_iter; fluid_model->SetTDState_T(val_temp, val_scalars); su2double Temperature = fluid_model->GetTemperature(); @@ -758,14 +832,14 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo su2double scalars[MAXNVAR], delta = 1e-3; for (auto iVar = 0u; iVar < nVar; iVar++) scalars[iVar] = scalar_solution[iVar]; bool outside = false; - scalars[I_PROGVAR] += delta; + scalars[idx_flame_front] += 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; - scalars[I_PROGVAR] += delta; + scalars[idx_flame_front] += delta; } - su2double pv_burnt = scalars[I_PROGVAR] - delta; + su2double pv_burnt = scalars[idx_flame_front] - 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; @@ -800,13 +874,13 @@ su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry* geometry, SU2_OMP_FOR_(schedule(static,omp_chunk_size) SU2_NOWAIT) for (auto iPoint = 0u; iPoint < nPointDomain; iPoint++) { - su2double pv_local = nodes->GetSolution(iPoint, I_PROGVAR); + su2double pv_local = nodes->GetSolution(iPoint, idx_progvar); su2double T_local = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint); /* Parallel projection of progress variable gradient against velocity */ su2double proj_grad_pv_u[MAXNDIM]={0}; for (auto iDim=0u; iDim < nDim; iDim++) { - su2double gradpv = nodes->GetGradient(iPoint, I_PROGVAR, iDim); + su2double gradpv = nodes->GetGradient(iPoint, idx_progvar, iDim); su2double val_u = flowNodes->GetVelocity(iPoint, iDim); proj_grad_pv_u[iDim] = gradpv * val_u * val_u / (flowNodes->GetVelocity2(iPoint) + EPS); } diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index 17e6e40807b5..2ce7ec50b157 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -238,6 +238,36 @@ void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; + /*--- Restart field of each variable, -1 for a variable that starts from its initial value. ---*/ + vector restart_field(nVar); + for (auto iVar = 0u; iVar < nVar; iVar++) restart_field[iVar] = skipVars + iVar; + + /*--- Flamelet scalars are found by name, a missing variance starts from its initial value. ---*/ + if (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET) { + const auto& options = config->GetFlameletParsedOptions(); + for (auto iVar = 0u; iVar < nVar; iVar++) { + const string& name = (iVar < options.n_control_vars) + ? options.controlling_variable_names[iVar] + : options.user_scalar_names[iVar - options.n_control_vars]; + const auto field = find_if(fields.begin(), fields.end(), [&name](string field_name) { + field_name.erase(remove(field_name.begin(), field_name.end(), '"'), field_name.end()); + return field_name == name; + }); + if (field != fields.end()) { + /*--- The fields start with the point index, which is not in the restart data. ---*/ + restart_field[iVar] = static_cast(distance(fields.begin(), field)) - 1; + } else if (static_cast(iVar) == options.idx_progvar_var || + static_cast(iVar) == options.idx_mixfrac_var) { + restart_field[iVar] = -1; + if (rank == MASTER_NODE) + cout << "The restart file has no field " << name << ", it starts from its initial value " + << config->GetSpecies_Init()[iVar] << "." << endl; + } else { + SU2_MPI::Error("The restart file " + restart_filename + " has no field " + name + ".", CURRENT_FUNCTION); + } + } + } + /*--- Load data from the restart into correct containers. ---*/ unsigned long counter = 0; @@ -251,9 +281,12 @@ void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi /*--- We need to store this point's data, so jump to the correct offset in the buffer of data from the restart file and load it. ---*/ - const auto index = counter * Restart_Vars[1] + skipVars; - for (auto iVar = 0u; iVar < nVar; iVar++) - nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]); + const auto index = counter * Restart_Vars[1]; + for (auto iVar = 0u; iVar < nVar; iVar++) { + const su2double value = (restart_field[iVar] >= 0) ? Restart_Data[index + restart_field[iVar]] + : config->GetSpecies_Init()[iVar]; + nodes->SetSolution(iPoint_Local, iVar, value); + } /*--- Increment the overall counter for how many points have been loaded. ---*/ counter++; diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index fe4214a965f2..37143747ef57 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -2062,7 +2062,7 @@ def main(): cfd_flamelet_ch4_cht.cfg_dir = "flamelet/03_laminar_premixed_ch4_flame_cht_cfd" cfd_flamelet_ch4_cht.cfg_file = "lam_prem_ch4_cht_cfd_master.cfg" cfd_flamelet_ch4_cht.test_iter = 5 - cfd_flamelet_ch4_cht.test_vals = [-7.223149, -6.338839, -8.334392, -10.410156, -3.171365, -10.884464, -12.615132, -6.614396] + cfd_flamelet_ch4_cht.test_vals = [-7.223149, -6.338839, -8.334392, -10.410156, -3.171365, -10.884424, -12.615132, -6.614396] cfd_flamelet_ch4_cht.timeout = 1600 cfd_flamelet_ch4_cht.multizone = True test_list.append(cfd_flamelet_ch4_cht) diff --git a/config_template.cfg b/config_template.cfg index 0ba2dc5b8924..e48b08c4c8a2 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1049,12 +1049,23 @@ LOOKUP_NAMES= (MolarWeightMix, Conductivity, HeatRelease, Diffusivity) % Names of the controlling variables used to evaluate the look-up table and for which % transported scalar variables are generated. These names should coincide with those -% listed in the look-up table. +% listed in the look-up table. Accepted names are ProgressVariable, EnthalpyTot, +% MixtureFraction, ProgressVariableVariance and MixtureFractionVariance. EnthalpyTot is +% required, together with ProgressVariable (premixed), MixtureFraction (non-premixed) or +% both. A variance requires its mean and KIND_TURB_MODEL= SST. The first two names are the +% table axes, the third one is the variable of the table levels. +% A variance is transported as such, but the table or network takes it normalized by its +% upper limit (mean - min)(max - mean), so with values in [0, 1]. For the mixture fraction +% min = 0 and max = 1. For the progress variable, min and max are the manifold variables +% ProgressVariableMin and ProgressVariableMax, which may depend on the other controlling +% variables but not on the variance. The variance is clipped to its upper limit. CONTROLLING_VARIABLE_NAMES= (ProgressVariable, EnthalpyTot) % Names of the source term variables as they are present in the look-up table. For % controlling variables without source terms, use "zero" or "null", similar to the -% option USER_SOURCE_NAMES. +% option USER_SOURCE_NAMES. For a variance, the variable is the mean of phi*S_phi, with +% S_phi the source of its mean phi; the chemical source of the variance is +% 2 (mean(phi*S_phi) - mean(phi) mean(S_phi)), or zero for NULL. CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL) % Method used to ignite the solution: @@ -1115,6 +1126,10 @@ SPARK_REACTION_RATES= (1000, 0, 0) % Approximate the thickness of the flame and suppress source terms and enhance diffusivity in cells with a lengthscale higher % than the flame thickness. THICKENED_FLAME_CORRECTION= YES +% +% Coefficient C_d of the dissipation C_d rho (eps/k) var of a transported scalar variance, +% with eps/k = beta* omega from the SST model. The production is 2 mu_t/Sc_t |grad(mean)|^2. +FLAMELET_VARIANCE_DISSIPATION_COEFF= 2.0 % Flame ignition temperature % Flame thickness calculations are performed when the maximum temperature in the flow solution exceeds this value.