diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index dffd3e85fdd..bc32b5407b1 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1224,7 +1224,13 @@ class CConfig { unsigned short nSpecies_Init; /*!< \brief Number of entries of SPECIES_INIT */ /*--- Additional flamelet solver options ---*/ - su2double flame_init[8]; /*!< \brief Initial solution parameters for flamelet solver.*/ + ///TODO: Add python wrapper initialization option + FLAMELET_INIT_TYPE flame_init_type = FLAMELET_INIT_TYPE::NONE; /*!< \brief Method for solution ignition for flamelet problems. */ + std::array flame_init; /*!< \brief Flame front initialization parameters. */ + std::array spark_init; /*!< \brief Spark ignition initialization parameters. */ + su2double* spark_reaction_rates; /*!< \brief Source terms for flamelet spark ignition option. */ + unsigned short nspark; /*!< \brief Number of source terms for spark initialization. */ + bool preferential_diffusion = false; /*!< \brief Preferential diffusion physics for flamelet solver.*/ /*--- lookup table ---*/ unsigned short n_scalars = 0; /*!< \brief Number of transported scalars for flamelet LUT approach. */ @@ -2138,13 +2144,44 @@ class CConfig { /*! * \brief Get the flame initialization. - * (x1,x2,x3) = flame offset. - * (x4,x5,x6) = flame normal, separating unburnt from burnt. + * (x1,x2,x3) = flame offset/spark center location. + * (x4,x5,x6) = flame normal, separating unburnt from burnt or + * spark radius, spark start iteration, spark duration. * (x7) = flame thickness, the length from unburnt to burnt conditions. * (x8) = flame burnt thickness, the length to stay at burnt conditions. - * \return Flame initialization for the flamelet model. + * \return Ignition initialization parameters for the flamelet model. + */ + const su2double* GetFlameInit() const { + switch (flame_init_type) + { + case FLAMELET_INIT_TYPE::FLAME_FRONT: + return flame_init.data(); + break; + case FLAMELET_INIT_TYPE::SPARK: + return spark_init.data(); + break; + default: + return nullptr; + break; + } + } + + /*! + * \brief Get species net reaction rates applied during spark ignition. + */ + const su2double* GetSpark() const { + return spark_reaction_rates; + } + + /*! + * \brief Preferential diffusion combustion problem. + */ + bool GetPreferentialDiffusion() const { return preferential_diffusion; } + + /*! + * \brief Define preferential diffusion combustion problem. */ - const su2double* GetFlameInit() const { return flame_init; } + inline void SetPreferentialDiffusion(bool input) { preferential_diffusion = input; } /*! * \brief Get the number of control variables for flamelet model. @@ -2190,6 +2227,11 @@ class CConfig { if (n_user_sources > 0) return user_source_names[i_user_source]; else return none; } + /*! + * \brief Get the ignition method used for combustion problems. + */ + FLAMELET_INIT_TYPE GetFlameletInitType() const { return flame_init_type; } + /*! * \brief Get the number of transported scalars for combustion. */ diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index 1c79535c5aa..230b009fc92 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/Common/include/containers/CLookUpTable.hpp @@ -389,7 +389,6 @@ class CLookUpTable { /*! * \brief Find the table levels with constant z-values directly above and below query val_z. * \param[in] val_CV3 - Value of controlling variable 3. - * \param[in] within_limits - Whether query point lies within table bounds. * \returns Pair of inclusion level indices (first = lower level index, second = upper level index). */ std::pair FindInclusionLevels(const su2double val_CV3); @@ -410,6 +409,11 @@ class CLookUpTable { return limits_table_x[i_level]; } + /*! + * \brief Check whether requested set of variables are included in the table. + */ + bool CheckForVariables(const std::vector& vars_to_check) const; + /*! * \brief Returns the index to the variable in the lookup table. * \param[in] nameVar - Variable name for which to retrieve the column index. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 7199a310193..919aefa65b9 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1333,11 +1333,38 @@ enum FLAMELET_SCALAR_SOURCES { * \brief Look-up operations for the flamelet scalar solver. */ enum FLAMELET_LOOKUP_OPS { - TD, /*!< \brief Thermochemical properties (temperature, density, diffusivity, etc.). */ + THERMO, /*!< \brief Thermochemical properties (temperature, density, diffusivity, etc.). */ + PREFDIF, /*!< \brief Preferential diffusion scalars. */ SOURCES, /*!< \brief Scalar source terms (controlling variables, passive species).*/ LOOKUP, /*!< \brief Passive look-up variables specified in config. */ }; +/*! + * \brief The preferential diffusion scalar indices for the preferential diffusion model. + */ +enum FLAMELET_PREF_DIFF_SCALARS { + I_BETA_PROGVAR, /*!< \brief Preferential diffusion scalar for the progress variable. */ + I_BETA_ENTH_THERMAL, /*!< \brief Preferential diffusion scalar for temperature. */ + I_BETA_ENTH, /*!< \brief Preferential diffusion scalar for total enthalpy. */ + I_BETA_MIXFRAC, /*!< \brief Preferential diffusion scalar for mixture fraction. */ + N_BETA_TERMS, /*!< \brief Total number of preferential diffusion scalars. */ +}; + +/*! + * \brief Flame initialization options for the flamelet solver. + */ +enum class FLAMELET_INIT_TYPE { + FLAME_FRONT, /*!< \brief Straight flame front. */ + SPARK, /*!< \brief Species reaction rate in a set location. */ + NONE, /*!< \brief No ignition, cold flow only. */ +}; + +static const MapType Flamelet_Init_Map = { + MakePair("NONE", FLAMELET_INIT_TYPE::NONE) + MakePair("FLAME_FRONT", FLAMELET_INIT_TYPE::FLAME_FRONT) + MakePair("SPARK", FLAMELET_INIT_TYPE::SPARK) +}; + /*! * \brief Types of subgrid scale models */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 18da76556c4..63810ab3e03 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -987,6 +987,7 @@ void CConfig::SetPointersNull() { Species_Init = nullptr; Species_Clipping_Min = nullptr; Species_Clipping_Max = nullptr; + spark_reaction_rates = nullptr; /*--- Moving mesh pointers ---*/ @@ -1376,14 +1377,23 @@ void CConfig::SetConfig_Options() { /*!\brief SPECIES_CLIPPING_MIN \n DESCRIPTION: Minimum values for scalar clipping \ingroup Config*/ addDoubleListOption("SPECIES_CLIPPING_MIN", nSpecies_Clipping_Min, Species_Clipping_Min); - /*!\brief FLAME_INIT \n DESCRIPTION: flame initialization using the flamelet model \ingroup Config*/ + /*!\brief FLAME_INIT_METHOD \n DESCRIPTION: Ignition method for flamelet solver \n DEFAULT: no ignition; cold flow only. */ + addEnumOption("FLAME_INIT_METHOD", flame_init_type, Flamelet_Init_Map, FLAMELET_INIT_TYPE::NONE); + /*!\brief FLAME_INIT \n DESCRIPTION: flame front initialization using the flamelet model \ingroup Config*/ /*--- flame offset (x,y,z) ---*/ flame_init[0] = 0.0; flame_init[1] = 0.0; flame_init[2] = 0.0; /*--- flame normal (nx, ny, nz) ---*/ flame_init[3] = 1.0; flame_init[4] = 0.0; flame_init[5] = 0.0; /*--- flame thickness (x) and flame burnt thickness (after this thickness, we have unburnt conditions again) ---*/ flame_init[6] = 0.5e-3; flame_init[7] = 1.0; - addDoubleArrayOption("FLAME_INIT", 8,flame_init); + addDoubleArrayOption("FLAME_INIT", 8,flame_init.begin()); + + /*!\brief SPARK_INIT \n DESCRIPTION: spark initialization using the flamelet model \ingroup Config*/ + for (auto iSpark=0u; iSpark<6; ++iSpark) spark_init[iSpark]=0; + addDoubleArrayOption("SPARK_INIT", 6, spark_init.begin()); + + /*!\brief SPARK_REACTION_RATES \n DESCRIPTION: Net source term values applied to species within spark area during spark ignition. \ingroup Config*/ + addDoubleListOption("SPARK_REACTION_RATES", nspark, spark_reaction_rates); /*--- Options related to mass diffusivity and thereby the species solver. ---*/ @@ -2136,6 +2146,9 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Names of the user scalar source terms. */ addStringListOption("USER_SOURCE_NAMES", n_user_sources, user_source_names); + /* DESCRIPTION: Enable preferential diffusion for FGM simulations. \n DEFAULT: false */ + addBoolOption("PREFERENTIAL_DIFFUSION", preferential_diffusion, false); + /*!\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*/ diff --git a/Common/src/containers/CLookUpTable.cpp b/Common/src/containers/CLookUpTable.cpp index 6efa79bcef9..16c042efba7 100644 --- a/Common/src/containers/CLookUpTable.cpp +++ b/Common/src/containers/CLookUpTable.cpp @@ -802,3 +802,13 @@ bool CLookUpTable::IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned l return (abs(area_tri - (area_0 + area_1 + area_2)) < area_tri * 1e-10); } + +bool CLookUpTable::CheckForVariables(const std::vector& vars_to_check) const { + for (const string& var_to_check : vars_to_check) { + if (!std::any_of(names_var.begin(), names_var.end(), + [var_to_check](const std::string& n) { return var_to_check == n; })) { + return false; + }; + } + return true; +} diff --git a/Docs/docmain.hpp b/Docs/docmain.hpp index 96743d0e371..5fefea50075 100644 --- a/Docs/docmain.hpp +++ b/Docs/docmain.hpp @@ -225,9 +225,3 @@ * \brief Classes for explicit (done by the programmer) vectorization (SIMD) of computations. * \ingroup Toolboxes */ - -/*! - * \defgroup Multi-Layer Perceptrons (MLP) - * \brief Data look up and interpolation via dense, feed-forward multi-layer perceptrons. - * \ingroup Toolboxes - */ diff --git a/SU2_CFD/include/fluid/CFluidFlamelet.hpp b/SU2_CFD/include/fluid/CFluidFlamelet.hpp index 12f370ad5eb..49273b2a9e6 100644 --- a/SU2_CFD/include/fluid/CFluidFlamelet.hpp +++ b/SU2_CFD/include/fluid/CFluidFlamelet.hpp @@ -42,9 +42,11 @@ class CFluidFlamelet final : public CFluidModel { enum LOOKUP_TD { TEMPERATURE, HEATCAPACITY, VISCOSITY, CONDUCTIVITY, DIFFUSIONCOEFFICIENT, MOLARWEIGHT, SIZE }; + su2double beta_progvar, beta_enth_thermal, beta_enth, beta_mixfrac; int rank; - bool include_mixture_fraction = false; /*!< \brief include mixture fraction in controlling variables. */ + bool include_mixture_fraction = false, /*!< \brief include mixture fraction in controlling variables. */ + 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. */ @@ -62,12 +64,15 @@ class CFluidFlamelet final : public CFluidModel { vector LUT_idx_TD, LUT_idx_Sources, - LUT_idx_LookUp; + LUT_idx_LookUp, + LUT_idx_PD; /*--- Class variables for the multi-layer perceptron method ---*/ #ifdef USE_MLPCPP + size_t n_betas; MLPToolbox::CLookUp_ANN* lookup_mlp; /*!< \brief Multi-layer perceptron collection. */ MLPToolbox::CIOMap* iomap_TD; /*!< \brief Input-output map for thermochemical properties. */ + 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; @@ -76,10 +81,14 @@ class CFluidFlamelet final : public CFluidModel { vector scalars_vector; vector varnames_TD, /*!< \brief Lookup names for thermodynamic state variables. */ - varnames_Sources, varnames_LookUp; + 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. */ vector val_vars_TD, /*!< \brief References to thermodynamic state variables. */ - val_vars_Sources, val_vars_LookUp; + 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. */ void PreprocessLookUp(CConfig* config); @@ -107,14 +116,14 @@ class CFluidFlamelet final : public CFluidModel { void SetTDState_T(su2double val_temperature, const su2double* val_scalars = nullptr) override; /*! - * \brief Evaluate the flamelet manifold. - * \param[in] input_scalar - scalar solution. - * \param[in] input_varnames - names of variables to evaluate. - * \param[in] output_refs - output data. - * \param[out] Extrapolation - scalar solution is within bounds (0) or out of bounds (1). + * \brief Evaluate data-set for flamelet simulations. + * \param[in] input_scalar - controlling variables used to interpolate manifold. + * \param[in] lookup_type - look-up operation to be performed (FLAMELET_LOOKUP_OPS) + * \param[in] output_refs - output variables where interpolated results are stored. + * \param[out] Extrapolation - query data is within manifold bounds (0) or out of bounds (1). */ - inline unsigned long EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, - vector& output_refs) override; + unsigned long EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, + vector& output_refs); /*! * \brief Check for out-of-bounds condition for data set interpolation. @@ -142,4 +151,10 @@ class CFluidFlamelet final : public CFluidModel { * \param[out] Mu - value of the laminar viscosity */ inline su2double GetLaminarViscosity() override { return Mu; } + + /*! + * \brief Preferential diffusion as relevant phenomenon in flamelet simulations. + * \return Inclusion of preferential diffusion model. + */ + inline bool GetPreferentialDiffusion() const override { return preferential_diffusion; } }; diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index 7fc767b17d2..29f0becdb0d 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/SU2_CFD/include/fluid/CFluidModel.hpp @@ -46,34 +46,34 @@ class CLookUpTable; */ class CFluidModel { protected: - su2double StaticEnergy{0.0}; /*!< \brief Internal Energy. */ - su2double Entropy{0.0}; /*!< \brief Entropy. */ - su2double Density{0.0}; /*!< \brief Density. */ - su2double Pressure{0.0}; /*!< \brief Pressure. */ - su2double SoundSpeed2{0.0}; /*!< \brief SpeedSound. */ - su2double Temperature{0.0}; /*!< \brief Temperature. */ - su2double dPdrho_e{0.0}; /*!< \brief DpDd_e. */ - su2double dPde_rho{0.0}; /*!< \brief DpDe_d. */ - su2double dTdrho_e{0.0}; /*!< \brief DTDd_e. */ - su2double dTde_rho{0.0}; /*!< \brief DTDe_d. */ - su2double dhdrho_P{0.0}; /*!< \brief DhDrho_p. */ - su2double dhdP_rho{0.0}; /*!< \brief DhDp_rho. */ - su2double dsdrho_P{0.0}; /*!< \brief DsDrho_p. */ - su2double dsdP_rho{0.0}; /*!< \brief DsDp_rho. */ - su2double Cp{0.0}; /*!< \brief Specific Heat Capacity at constant pressure. */ - su2double Cv{0.0}; /*!< \brief Specific Heat Capacity at constant volume. */ - su2double Mu{0.0}; /*!< \brief Laminar viscosity. */ - su2double Mu_Turb{0.0}; /*!< \brief Eddy viscosity provided by a turbulence model. */ - su2double dmudrho_T{0.0}; /*!< \brief Partial derivative of viscosity w.r.t. density. */ - su2double dmudT_rho{0.0}; /*!< \brief Partial derivative of viscosity w.r.t. temperature. */ - su2double Kt{0.0}; /*!< \brief Thermal conductivity. */ - su2double dktdrho_T{0.0}; /*!< \brief Partial derivative of conductivity w.r.t. density. */ - su2double dktdT_rho{0.0}; /*!< \brief Partial derivative of conductivity w.r.t. temperature. */ - su2double mass_diffusivity{0.0}; /*!< \brief Mass Diffusivity */ + su2double StaticEnergy{0.0}; /*!< \brief Internal Energy. */ + su2double Entropy{0.0}; /*!< \brief Entropy. */ + su2double Density{0.0}; /*!< \brief Density. */ + su2double Pressure{0.0}; /*!< \brief Pressure. */ + su2double SoundSpeed2{0.0}; /*!< \brief SpeedSound. */ + su2double Temperature{0.0}; /*!< \brief Temperature. */ + su2double dPdrho_e{0.0}; /*!< \brief DpDd_e. */ + su2double dPde_rho{0.0}; /*!< \brief DpDe_d. */ + su2double dTdrho_e{0.0}; /*!< \brief DTDd_e. */ + su2double dTde_rho{0.0}; /*!< \brief DTDe_d. */ + su2double dhdrho_P{0.0}; /*!< \brief DhDrho_p. */ + su2double dhdP_rho{0.0}; /*!< \brief DhDp_rho. */ + su2double dsdrho_P{0.0}; /*!< \brief DsDrho_p. */ + su2double dsdP_rho{0.0}; /*!< \brief DsDp_rho. */ + su2double Cp{0.0}; /*!< \brief Specific Heat Capacity at constant pressure. */ + su2double Cv{0.0}; /*!< \brief Specific Heat Capacity at constant volume. */ + su2double Mu{0.0}; /*!< \brief Laminar viscosity. */ + su2double Mu_Turb{0.0}; /*!< \brief Eddy viscosity provided by a turbulence model. */ + su2double dmudrho_T{0.0}; /*!< \brief Partial derivative of viscosity w.r.t. density. */ + su2double dmudT_rho{0.0}; /*!< \brief Partial derivative of viscosity w.r.t. temperature. */ + su2double Kt{0.0}; /*!< \brief Thermal conductivity. */ + su2double dktdrho_T{0.0}; /*!< \brief Partial derivative of conductivity w.r.t. density. */ + su2double dktdT_rho{0.0}; /*!< \brief Partial derivative of conductivity w.r.t. temperature. */ + su2double mass_diffusivity{0.0}; /*!< \brief Mass Diffusivity */ unique_ptr LaminarViscosity; /*!< \brief Laminar Viscosity Model */ unique_ptr ThermalConductivity; /*!< \brief Thermal Conductivity Model */ - unique_ptr MassDiffusivity; /*!< \brief Mass Diffusivity Model */ + unique_ptr MassDiffusivity; /*!< \brief Mass Diffusivity Model */ /*! * \brief Instantiate the right type of viscosity model based on config. @@ -144,10 +144,16 @@ class CFluidModel { virtual inline unsigned short GetNScalars() const { return 0; } /*! - * \brief Evaluate data manifold for flamelet or data-driven fluid problems. - * \param[in] input - input data for manifold regression. + * \brief Evaluate data-set for flamelet or data-driven fluid simulations. + * \param[in] input_scalar - data manifold query data. + * \param[in] lookup_type - look-up operation to be performed. + * \param[in] output_refs - output variables where interpolated results are stored. + * \param[out] Extrapolation - query data is within manifold bounds (0) or out of bounds (1). */ - virtual unsigned long EvaluateDataSet(const vector &input_scalar, unsigned short lookup_type, vector &output_refs) { return 0; } + virtual unsigned long EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, + vector& output_refs) { + return 0; + } /*! * \brief Get fluid dynamic viscosity. @@ -331,7 +337,7 @@ class CFluidModel { * \brief Virtual member. * \param[in] T - Temperature value at the point. */ - virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) { } + virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) {} /*! * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity. @@ -344,6 +350,12 @@ class CFluidModel { */ virtual unsigned long GetExtrapolation() const { return 0; } + /*! + * \brief Get the state of the Preferential diffusion model for flamelet simulations. + * \return True if preferential diffusion model is active, false otherwise. + */ + virtual bool GetPreferentialDiffusion() const { return false; } + /*! * \brief Get number of Newton solver iterations. * \return Newton solver iteration count at termination. diff --git a/SU2_CFD/include/solvers/CHeatSolver.hpp b/SU2_CFD/include/solvers/CHeatSolver.hpp index e82c96881d7..627983ec096 100644 --- a/SU2_CFD/include/solvers/CHeatSolver.hpp +++ b/SU2_CFD/include/solvers/CHeatSolver.hpp @@ -103,8 +103,8 @@ class CHeatSolver final : public CScalarSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - inline void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) override { + inline void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override { const CVariable* flow_nodes = flow ? solver_container[FLOW_SOL]->GetNodes() : nullptr; const su2double const_diffusivity = config->GetThermalDiffusivity(); diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index ca84ea7c644..ca4ce3e2483 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -99,9 +99,9 @@ class CScalarSolver : public CSolver { * \param[in] config - Definition of the particular problem. */ template - FORCEINLINE void Viscous_Residual_impl(const SolverSpecificNumericsFunc& SolverSpecificNumerics, unsigned long iEdge, - CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, - CConfig* config) { + FORCEINLINE void Viscous_Residual_impl(const SolverSpecificNumericsFunc& SolverSpecificNumerics, const unsigned long iEdge, + const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, + const CConfig* config) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); CFlowVariable* flowNodes = solver_container[FLOW_SOL] ? su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()) : nullptr; @@ -307,8 +307,8 @@ class CScalarSolver : public CSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - inline virtual void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { + inline virtual void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { /*--- Define an empty object for solver specific numerics contribution. In case there are none, this default *--- implementation will be called ---*/ auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) {}; diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index e821e1d6a9e..1d1bd498230 100644 --- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp @@ -37,7 +37,6 @@ */ class CSpeciesFlameletSolver final : public CSpeciesSolver { private: - vector conjugate_var; /*!< \brief CHT variables for each boundary and vertex. */ bool include_mixture_fraction = false; /*!< \brief include mixture fraction as a controlling variable. */ /*! * \brief Compute the preconditioner for low-Mach flows. @@ -63,8 +62,8 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { * \param[in] val_enth_out - pointer to output enthalpy variable. * \param[out] Converged - 0 if Newton solver converged, 1 if not. */ - unsigned long GetEnthFromTemp(CFluidModel* fluid_model, su2double const val_temp, - const su2double* scalar_solution, su2double* val_enth_out); + unsigned long GetEnthFromTemp(CFluidModel* fluid_model, su2double const val_temp, const su2double* scalar_solution, + su2double* val_enth_out); /*! * \brief Find maximum progress variable value within the manifold for the current solution. @@ -97,13 +96,23 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { unsigned long SetScalarLookUps(const CConfig* config, CFluidModel* fluid_model_local, unsigned long iPoint, const vector& scalars); + /*! + * \brief Retrieve the preferential diffusion scalar values from manifold. + * \param[in] config - definition of particular problem. + * \param[in] fluid_model_local - pointer to flamelet fluid model. + * \param[in] iPoint - node ID. + * \param[in] scalars - local scalar solution. + * \return - within manifold bounds (0) or outside manifold bounds (1). + */ + unsigned long SetPreferentialDiffusionScalars(const CConfig* config, CFluidModel* fluid_model_local, + unsigned long iPoint, const vector& scalars); + public: /*! - * \brief Constructor. + * \brief Define a Flamelet Generated Manifold species solver. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. * \param[in] iMesh - Index of the mesh in multigrid computations. - * \param[in] FluidModel */ CSpeciesFlameletSolver(CGeometry* geometry, CConfig* config, unsigned short iMesh); @@ -175,4 +184,16 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { */ void BC_ConjugateHeat_Interface(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, CConfig* config, unsigned short val_marker) override; + + /*! + * \brief Compute the fluxes due to viscous and preferential diffusion effects of the flamelet species at a particular edge. + * \param[in] iEdge - Edge for which we want to compute the flux + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \note Calls a generic implementation after defining a SolverSpecificNumerics object. + */ + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, + const CConfig* config) final; }; diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 1ee78d90286..4fbb9419117 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -89,8 +89,8 @@ class CSpeciesSolver : public CScalarSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, - CConfig* config) final; + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, + const CConfig* config) override; /*! * \brief Impose the inlet boundary condition. diff --git a/SU2_CFD/include/solvers/CTransLMSolver.hpp b/SU2_CFD/include/solvers/CTransLMSolver.hpp index 220f2b01ed7..98595abb266 100644 --- a/SU2_CFD/include/solvers/CTransLMSolver.hpp +++ b/SU2_CFD/include/solvers/CTransLMSolver.hpp @@ -98,8 +98,8 @@ class CTransLMSolver final : public CTurbSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) override; + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override; /*! * \brief Source term computation. diff --git a/SU2_CFD/include/solvers/CTurbSASolver.hpp b/SU2_CFD/include/solvers/CTurbSASolver.hpp index cb0307817e7..c79649799cc 100644 --- a/SU2_CFD/include/solvers/CTurbSASolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSASolver.hpp @@ -125,8 +125,8 @@ class CTurbSASolver final : public CTurbSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) override; + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override; /*! * \brief Source term computation. diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp index f4cf89d38d0..bdb4227a191 100644 --- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp @@ -107,8 +107,8 @@ class CTurbSSTSolver final : public CTurbSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) override; + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override; /*! * \brief Source term computation. diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index 3e7c67c2c72..18dbcd2a43a 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -25,6 +25,8 @@ * License along with SU2. If not, see . */ +#include +#include #include "../include/fluid/CFluidFlamelet.hpp" #include "../../../Common/include/containers/CLookUpTable.hpp" #if defined(HAVE_MLPCPP) @@ -47,7 +49,6 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati cout << "Number of user scalars: " << n_user_scalars << endl; cout << "Number of control variables: " << n_control_vars << endl; } - scalars_vector.resize(n_scalars); table_scalar_names.resize(n_scalars); @@ -93,31 +94,31 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati Pressure = value_pressure_operating; PreprocessLookUp(config); + + if (rank == MASTER_NODE) { + cout << "Preferential diffusion: " << (preferential_diffusion ? "Enabled" : "Disabled") << endl; + } } CFluidFlamelet::~CFluidFlamelet() { - switch (Kind_DataDriven_Method) { - case ENUM_DATADRIVEN_METHOD::LUT: - delete look_up_table; - break; - case ENUM_DATADRIVEN_METHOD::MLP: + if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT) + delete look_up_table; #ifdef USE_MLPCPP - delete iomap_TD; - delete iomap_Sources; - delete iomap_LookUp; - delete lookup_mlp; + if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::MLP) { + delete iomap_TD; + delete iomap_Sources; + delete iomap_LookUp; + delete lookup_mlp; + if (preferential_diffusion) delete iomap_PD; + } #endif - break; - default: - break; - } } void CFluidFlamelet::SetTDState_T(su2double val_temperature, const su2double* val_scalars) { for (auto iVar = 0u; iVar < n_scalars; iVar++) scalars_vector[iVar] = val_scalars[iVar]; /*--- Add all quantities and their names to the look up vectors. ---*/ - EvaluateDataSet(scalars_vector, FLAMELET_LOOKUP_OPS::TD, val_vars_TD); + EvaluateDataSet(scalars_vector, FLAMELET_LOOKUP_OPS::THERMO, val_vars_TD); Temperature = val_vars_TD[LOOKUP_TD::TEMPERATURE]; Cp = val_vars_TD[LOOKUP_TD::HEATCAPACITY]; @@ -183,9 +184,51 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { /*--- Passive look-up terms ---*/ size_t n_lookups = config->GetNLookups(); - varnames_LookUp.resize(n_lookups); - val_vars_LookUp.resize(n_lookups); - for (auto iLookup = 0u; iLookup < n_lookups; iLookup++) varnames_LookUp[iLookup] = config->GetLookupName(iLookup); + if (n_lookups == 0) { + varnames_LookUp.resize(1); + val_vars_LookUp.resize(1); + varnames_LookUp[0] = "NULL"; + } else { + varnames_LookUp.resize(n_lookups); + val_vars_LookUp.resize(n_lookups); + for (auto iLookup = 0u; iLookup < n_lookups; iLookup++) varnames_LookUp[iLookup] = config->GetLookupName(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); + + 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"; + + 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; + + preferential_diffusion = config->GetPreferentialDiffusion(); + switch (Kind_DataDriven_Method) { + case ENUM_DATADRIVEN_METHOD::LUT: + preferential_diffusion = look_up_table->CheckForVariables(varnames_PD); + break; + case ENUM_DATADRIVEN_METHOD::MLP: +#ifdef USE_MLPCPP + n_betas = 0; + for (auto iMLP = 0u; iMLP < config->GetNDataDriven_Files(); iMLP++) { + auto outputMap = lookup_mlp->FindVariableIndices(iMLP, varnames_PD, false); + n_betas += outputMap.size(); + } + preferential_diffusion = (n_betas == varnames_PD.size()); +#endif + break; + default: + break; + } + + if (!preferential_diffusion && config->GetPreferentialDiffusion()) + SU2_MPI::Error("Preferential diffusion scalars not included in flamelet manifold.", CURRENT_FUNCTION); if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::MLP) { #ifdef USE_MLPCPP @@ -195,6 +238,10 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { lookup_mlp->PairVariableswithMLPs(*iomap_TD); lookup_mlp->PairVariableswithMLPs(*iomap_Sources); lookup_mlp->PairVariableswithMLPs(*iomap_LookUp); + if (preferential_diffusion) { + iomap_PD = new MLPToolbox::CIOMap(controlling_variable_names, varnames_PD); + lookup_mlp->PairVariableswithMLPs(*iomap_PD); + } #endif } else { for (auto iVar=0u; iVar < varnames_TD.size(); iVar++) { @@ -210,37 +257,52 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { LUT_idx_Sources.push_back(LUT_idx); } for (auto iVar=0u; iVar < varnames_LookUp.size(); iVar++) { - LUT_idx_LookUp.push_back(look_up_table->GetIndexOfVar(varnames_LookUp[iVar])); + unsigned long LUT_idx; + if (noSource(varnames_LookUp[iVar])) + LUT_idx = look_up_table->GetNullIndex(); + else + LUT_idx = look_up_table->GetIndexOfVar(varnames_LookUp[iVar]); + LUT_idx_LookUp.push_back(LUT_idx); + } + if (preferential_diffusion) { + for (auto iVar=0u; iVar < varnames_PD.size(); iVar++) { + LUT_idx_PD.push_back(look_up_table->GetIndexOfVar(varnames_PD[iVar])); + } } } } 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 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 varnames; vector val_vars; vector refs_vars; vector LUT_idx; switch (lookup_type) { - case FLAMELET_LOOKUP_OPS::TD: - varnames = varnames_TD; + case FLAMELET_LOOKUP_OPS::THERMO: LUT_idx = LUT_idx_TD; #ifdef USE_MLPCPP iomap_Current = iomap_TD; +#endif + break; + case FLAMELET_LOOKUP_OPS::PREFDIF: + LUT_idx = LUT_idx_PD; +#ifdef USE_MLPCPP + iomap_Current = iomap_PD; #endif break; case FLAMELET_LOOKUP_OPS::SOURCES: - varnames = varnames_Sources; LUT_idx = LUT_idx_Sources; #ifdef USE_MLPCPP iomap_Current = iomap_Sources; #endif break; case FLAMELET_LOOKUP_OPS::LOOKUP: - varnames = varnames_LookUp; LUT_idx = LUT_idx_LookUp; #ifdef USE_MLPCPP iomap_Current = iomap_LookUp; @@ -249,13 +311,14 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca default: break; } - if (output_refs.size() != varnames.size()) - SU2_MPI::Error(string("Output vector size incompatible with manifold lookup operation."), CURRENT_FUNCTION); + /*--- Add all quantities and their names to the look up vectors. ---*/ bool inside; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: + 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); } else { @@ -274,6 +337,7 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca default: break; } - + 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/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 66b90d3e39b..661c9543024 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -68,9 +68,19 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver unsigned short RunTime_EqSystem, bool Output) { unsigned long n_not_in_domain_local = 0, n_not_in_domain_global = 0; vector scalars_vector(nVar); - + unsigned long spark_iter_start, spark_duration; + bool ignition = false; auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + /*--- Retrieve spark ignition parameters for spark-type ignition. ---*/ + if ((config->GetFlameletInitType() == FLAMELET_INIT_TYPE::SPARK) && !config->GetRestart()) { + auto spark_init = config->GetFlameInit(); + spark_iter_start = ceil(spark_init[4]); + spark_duration = ceil(spark_init[5]); + unsigned long iter = config->GetMultizone_Problem() ? config->GetOuterIter() : config->GetInnerIter(); + ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)) && !config->GetRestart()); + } + SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);) SU2_OMP_FOR_STAT(omp_chunk_size) @@ -81,19 +91,35 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver /*--- Compute total source terms from the production and consumption. ---*/ unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector); + + if (ignition) { + /*--- Apply source terms within spark radius. ---*/ + su2double dist_from_center = 0, + spark_radius = config->GetFlameInit()[3]; + dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), config->GetFlameInit()); + if (dist_from_center < pow(spark_radius,2)) { + for (auto iVar = 0u; iVar < nVar; iVar++) + nodes->SetScalarSource(i_point, iVar, nodes->GetScalarSources(i_point)[iVar] + config->GetSpark()[iVar]); + } + } + nodes->SetTableMisses(i_point, misses); n_not_in_domain_local += misses; /*--- Obtain passive look-up scalars. ---*/ SetScalarLookUps(config, fluid_model_local, i_point, scalars_vector); /*--- Set mass diffusivity based on thermodynamic state. ---*/ - su2double T = flowNodes->GetTemperature(i_point); + auto T = flowNodes->GetTemperature(i_point); fluid_model_local->SetTDState_T(T, scalars); /*--- set the diffusivity in the fluid model to the diffusivity obtained from the lookup table ---*/ for (auto i_scalar = 0u; i_scalar < nVar; ++i_scalar) { nodes->SetDiffusivity(i_point, fluid_model_local->GetMassDiffusivity(i_scalar), i_scalar); } + /*--- Obtain preferential diffusion scalar values. ---*/ + if (config->GetPreferentialDiffusion()) + SetPreferentialDiffusionScalars(config, fluid_model_local, i_point, scalars_vector); + if (!Output) LinSysRes.SetBlock_Zero(i_point); } END_SU2_OMP_FOR @@ -102,6 +128,20 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver SU2_MPI::GetComm()); 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 (config->GetPreferentialDiffusion()) { + 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; + } + } /*--- Clear Residual and Jacobian. Upwind second order reconstruction and gradients ---*/ CommonPreprocessing(geometry, config, Output); } @@ -115,13 +155,22 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** cout << "Initializing progress variable and total enthalpy (using temperature)" << endl; } - const su2double* flame_init = config->GetFlameInit(); - const su2double flame_offset[3] = {flame_init[0], flame_init[1], flame_init[2]}; - const su2double flame_normal[3] = {flame_init[3], flame_init[4], flame_init[5]}; - const su2double flame_thickness = flame_init[6]; - const su2double flame_burnt_thickness = flame_init[7]; + su2double flame_offset[3] = {0, 0, 0}, flame_normal[3] = {0, 0, 0}, flame_thickness = 0, flame_burnt_thickness = 0, + flamenorm = 0; + bool flame_front_ignition = (config->GetFlameletInitType() == FLAMELET_INIT_TYPE::FLAME_FRONT); + + if (flame_front_ignition) { + /*--- Collect flame front ignition parameters. ---*/ + auto flame_init = config->GetFlameInit(); + for (auto iDim = 0u; iDim < 3; ++iDim) { + flame_offset[iDim] = flame_init[iDim]; + flame_normal[iDim] = flame_init[3 + iDim]; + } + flame_thickness = flame_init[6]; + flame_burnt_thickness = flame_init[7]; + flamenorm = GeometryToolbox::Norm(nDim, flame_normal); + } - const su2double flamenorm = GeometryToolbox::Norm(nDim, flame_normal); const su2double temp_inlet = config->GetInc_Temperature_Init(); su2double enth_inlet = config->GetSpecies_Init()[I_ENTH]; @@ -134,6 +183,19 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** const auto& cv_name = config->GetControllingVariableName(iCV); cout << "initial condition: " << cv_name << " = " << config->GetSpecies_Init()[iCV] << endl; } + switch (config->GetFlameletInitType()) { + case FLAMELET_INIT_TYPE::FLAME_FRONT: + cout << "Ignition with a straight flame front" << endl; + break; + case FLAMELET_INIT_TYPE::SPARK: + cout << "Ignition with an artificial spark" << endl; + break; + case FLAMELET_INIT_TYPE::NONE: + cout << "No solution ignition (cold flow)" << endl; + break; + default: + break; + } } CFluidModel* fluid_model_local; @@ -157,37 +219,40 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** for (unsigned long i_point = 0; i_point < nPoint; i_point++) { auto coords = geometry[i_mesh]->nodes->GetCoord(i_point); - /*--- Determine if point is above or below the plane, assuming the normal - is pointing towards the burned region. ---*/ - point_loc = 0.0; - for (unsigned short i_dim = 0; i_dim < geometry[i_mesh]->GetnDim(); i_dim++) { - point_loc += flame_normal[i_dim] * (coords[i_dim] - flame_offset[i_dim]); - } - - /*--- Compute the exact distance from point to plane. ---*/ - point_loc = point_loc / flamenorm; - - /* --- Unburnt region upstream of the flame. --- */ - if (point_loc <= 0) { - scalar_init[I_PROGVAR] = 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; - 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; - n_points_burnt_local++; - - /* --- Unburnt region downstream of the flame and burnt region. --- */ + if (flame_front_ignition) { + /*--- Determine if point is above or below the plane, assuming the normal + is pointing towards the burned region. ---*/ + point_loc = 0.0; + for (unsigned short i_dim = 0; i_dim < geometry[i_mesh]->GetnDim(); i_dim++) { + point_loc += flame_normal[i_dim] * (coords[i_dim] - flame_offset[i_dim]); + } + + /*--- Compute the exact distance from point to plane. ---*/ + point_loc = point_loc / flamenorm; + + /* --- Unburnt region upstream of the flame. --- */ + if (point_loc <= 0) { + scalar_init[I_PROGVAR] = 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; + 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; + n_points_burnt_local++; + + /* --- Unburnt region downstream of the flame and burnt region. --- */ + } else { + scalar_init[I_PROGVAR] = prog_unburnt; + n_points_unburnt_local++; + } } else { scalar_init[I_PROGVAR] = prog_unburnt; - n_points_unburnt_local++; } - /* --- Perform manifold evaluation to check whether initial scalar is out-of-bounds. --- */ fluid_model_local->SetTDState_T(temp_inlet, scalar_init); n_not_in_domain_local += fluid_model_local->GetExtrapolation(); @@ -225,9 +290,11 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** if (rank == MASTER_NODE) { cout << endl; - cout << " Number of points in unburnt region: " << n_points_unburnt_global << "." << endl; - cout << " Number of points in burnt region : " << n_points_burnt_global << "." << endl; - cout << " Number of points in flame zone : " << n_points_flame_global << "." << endl; + if (flame_front_ignition) { + cout << " Number of points in unburnt region: " << n_points_unburnt_global << "." << endl; + cout << " Number of points in burnt region : " << n_points_burnt_global << "." << endl; + cout << " Number of points in flame zone : " << n_points_flame_global << "." << endl; + } if (n_not_in_domain_global > 0) cout << " Initial condition: Number of points outside of table domain: " << n_not_in_domain_global << " !!!" @@ -331,8 +398,8 @@ void CSpeciesFlameletSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_cont 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; - END_SU2_OMP_FOR } + END_SU2_OMP_FOR /*--- Call the general inlet boundary condition implementation. ---*/ CSpeciesSolver::BC_Inlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); @@ -343,10 +410,8 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo CConfig* config, unsigned short val_marker, bool cht_mode) { const bool implicit = config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT; const string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - su2double temp_wall; CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); - su2double enth_wall; unsigned long n_not_iterated = 0; /*--- Loop over all the vertices on this boundary marker. ---*/ @@ -354,6 +419,7 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo for (unsigned long iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { unsigned long iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + su2double temp_wall, enth_wall; if (cht_mode) temp_wall = solver_container[FLOW_SOL]->GetConjugateHeatVariable(val_marker, iVertex, 0); else @@ -402,7 +468,7 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo su2double dist_ij = sqrt(dist_ij_2); /*--- Compute the normal gradient in temperature using Twall. ---*/ - + ///TODO: Account for preferential diffusion in computation of the heat flux su2double dTdn = -(flowNodes->GetTemperature(Point_Normal) - temp_wall) / dist_ij; /*--- Get thermal conductivity. ---*/ @@ -442,13 +508,13 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF vector table_sources(config->GetNControlVars() + 2 * config->GetNUserScalars()); 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 makes sure. --- */ vector source_scalar(config->GetNScalars()); for (auto iCV = 0u; iCV < config->GetNControlVars(); iCV++) source_scalar[iCV] = table_sources[iCV]; - source_scalar[I_PROGVAR] = fmax(EPS, source_scalar[I_PROGVAR]); /*--- Source term for the auxiliary species transport equations. ---*/ for (size_t i_aux = 0; i_aux < config->GetNUserScalars(); i_aux++) { @@ -466,17 +532,212 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF unsigned long CSpeciesFlameletSolver::SetScalarLookUps(const CConfig* config, CFluidModel* fluid_model_local, unsigned long iPoint, const vector& scalars) { - /*--- Compute total source terms from the production and consumption. ---*/ + /*--- Retrieve the passive look-up variables from the manifold. ---*/ + unsigned long misses{0}; + /*--- Skip if no passive look-ups are listed ---*/ + if (config->GetNLookups() > 0) { + vector lookup_scalar(config->GetNLookups()); + misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::LOOKUP, lookup_scalar); + + for (auto i_lookup = 0u; i_lookup < config->GetNLookups(); i_lookup++) { + nodes->SetLookupScalar(iPoint, lookup_scalar[i_lookup], i_lookup); + } + } + + return misses; +} - vector lookup_scalar(config->GetNLookups()); - unsigned long misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::LOOKUP, lookup_scalar); +unsigned long CSpeciesFlameletSolver::SetPreferentialDiffusionScalars(const CConfig* config, + CFluidModel* fluid_model_local, + unsigned long iPoint, + const vector& scalars) { + /*--- Retrieve the preferential diffusion scalar values from the manifold. ---*/ - for (auto i_lookup = 0u; i_lookup < config->GetNLookups(); i_lookup++) { - nodes->SetLookupScalar(iPoint, lookup_scalar[i_lookup], i_lookup); + vector beta_scalar(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); + unsigned long misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::PREFDIF, beta_scalar); + + 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]); } return misses; } +void CSpeciesFlameletSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { + /*--- Overloaded viscous residual method which accounts for preferential diffusion. ---*/ + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT), + PreferentialDiffusion = config->GetPreferentialDiffusion(); + + /*--- Points in edge ---*/ + auto iPoint = geometry->edges->GetNode(iEdge, 0); + auto jPoint = geometry->edges->GetNode(iEdge, 1); + + auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) { + /*--- Mass diffusivity coefficients. ---*/ + + numerics->SetDiffusionCoeff(nodes->GetDiffusivity(iPoint), nodes->GetDiffusivity(jPoint)); + }; + + /*--- Regular viscous scalar residual computation. ---*/ + + Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config); + + /*--- Viscous residual due to preferential diffusion ---*/ + if (PreferentialDiffusion) { + CFlowVariable* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + + su2double scalar_i[MAXNVAR] = {0}, + scalar_j[MAXNVAR] = {0}, + diff_coeff_beta_i[MAXNVAR] = {0}, + diff_coeff_beta_j[MAXNVAR] = {0}; + + // Number of active transport scalars + const auto n_CV = config->GetNControlVars(); + + su2activematrix scalar_grad_i(MAXNVAR, MAXNDIM), scalar_grad_j(MAXNVAR, MAXNDIM); + /*--- Looping over spatial dimensions to fill in the diffusion scalar gradients. ---*/ + /*--- The scalar gradient is subtracted to account for regular viscous diffusion. ---*/ + for (auto iScalar = 0u; iScalar < n_CV; ++iScalar) { + for (auto iDim = 0u; iDim < nDim; ++iDim) { + switch (iScalar) { + case I_PROGVAR: + scalar_grad_i[iScalar][iDim] = + nodes->GetAuxVarGradient(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR, iDim) - + nodes->GetGradient(iPoint, I_PROGVAR, iDim); + scalar_grad_j[iScalar][iDim] = + nodes->GetAuxVarGradient(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR, iDim) - + nodes->GetGradient(jPoint, I_PROGVAR, iDim); + break; + case I_ENTH: + scalar_grad_i[iScalar][iDim] = + nodes->GetAuxVarGradient(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH, iDim) - + nodes->GetGradient(iPoint, I_ENTH, iDim); + scalar_grad_j[iScalar][iDim] = + nodes->GetAuxVarGradient(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH, iDim) - + nodes->GetGradient(jPoint, I_ENTH, iDim); + break; + case I_MIXFRAC: + scalar_grad_i[iScalar][iDim] = + nodes->GetAuxVarGradient(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC, iDim) - + nodes->GetGradient(iPoint, I_MIXFRAC, iDim); + scalar_grad_j[iScalar][iDim] = + nodes->GetAuxVarGradient(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC, iDim) - + nodes->GetGradient(jPoint, I_MIXFRAC, iDim); + break; + default: + break; + } + } + } + /*--- No preferential diffusion modification for passive species. ---*/ + for (auto iScalar = n_CV; iScalar < nVar; ++iScalar) { + for (auto iDim = 0u; iDim < nDim; ++iDim) { + scalar_grad_i[iScalar][iDim] = 0; + scalar_grad_j[iScalar][iDim] = 0; + } + } + + for (auto iScalar = 0u; iScalar < n_CV; ++iScalar) { + /*--- Filling in the preferential diffusion scalars (beta_pv, beta_h2, beta_Z). ---*/ + switch (iScalar) { + case I_PROGVAR: + scalar_i[iScalar] = nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR) - + nodes->GetSolution(iPoint, iScalar); + scalar_j[iScalar] = nodes->GetAuxVar(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR) - + nodes->GetSolution(jPoint, iScalar); + break; + case I_ENTH: + scalar_i[iScalar] = + nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH) - nodes->GetSolution(iPoint, iScalar); + scalar_j[iScalar] = + nodes->GetAuxVar(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH) - nodes->GetSolution(jPoint, iScalar); + break; + case I_MIXFRAC: + scalar_i[iScalar] = nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC) - + nodes->GetSolution(iPoint, iScalar); + scalar_j[iScalar] = nodes->GetAuxVar(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC) - + nodes->GetSolution(jPoint, iScalar); + break; + default: + break; + } + diff_coeff_beta_i[iScalar] = nodes->GetDiffusivity(iPoint, iScalar); + diff_coeff_beta_j[iScalar] = nodes->GetDiffusivity(jPoint, iScalar); + } + + for (auto iScalar = n_CV; iScalar < nVar; ++iScalar) { + scalar_i[iScalar] = 0; + scalar_j[iScalar] = 0; + diff_coeff_beta_i[iScalar] = 0; + diff_coeff_beta_j[iScalar] = 0; + } + + numerics->SetScalarVar(scalar_i, scalar_j); + + numerics->SetScalarVarGradient(CMatrixView(scalar_grad_i), CMatrixView(scalar_grad_j)); + + numerics->SetDiffusionCoeff(diff_coeff_beta_i, diff_coeff_beta_j); + + /*--- Computing first preferential residual component. ---*/ + auto residual_PD = numerics->ComputeResidual(config); + + if (ReducerStrategy) { + EdgeFluxes.SubtractBlock(iEdge, residual_PD); + + if (implicit) Jacobian.UpdateBlocksSub(iEdge, residual_PD.jacobian_i, residual_PD.jacobian_j); + } else { + LinSysRes.SubtractBlock(iPoint, residual_PD); + LinSysRes.AddBlock(jPoint, residual_PD); + /*--- Set implicit computation ---*/ + if (implicit) Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual_PD.jacobian_i, residual_PD.jacobian_j); + } + + /* Computing the second preferential diffusion terms due to heat flux */ + for (auto iScalar = 0u; iScalar < nVar; ++iScalar) { + for (auto iDim = 0u; iDim < nDim; ++iDim) { + if (iScalar == I_ENTH) { + /* Setting the temperature gradient */ + scalar_grad_i[iScalar][iDim] = flowNodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim); + scalar_grad_j[iScalar][iDim] = flowNodes->GetGradient_Primitive(jPoint, prim_idx.Temperature(), iDim); + } else { + scalar_grad_i[iScalar][iDim] = 0; + scalar_grad_j[iScalar][iDim] = 0; + } + } + + if (iScalar == I_ENTH) { + scalar_i[iScalar] = flowNodes->GetTemperature(iPoint); + scalar_j[iScalar] = flowNodes->GetTemperature(jPoint); + diff_coeff_beta_i[iScalar] = nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL) * + nodes->GetDiffusivity(iPoint, iScalar); + diff_coeff_beta_j[iScalar] = nodes->GetAuxVar(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL) * + nodes->GetDiffusivity(jPoint, iScalar); + } else { + scalar_i[iScalar] = 0; + scalar_j[iScalar] = 0; + diff_coeff_beta_i[iScalar] = 0; + diff_coeff_beta_j[iScalar] = 0; + } + } + + numerics->SetScalarVar(scalar_i, scalar_j); + + numerics->SetScalarVarGradient(CMatrixView(scalar_grad_i), CMatrixView(scalar_grad_j)); + + numerics->SetDiffusionCoeff(diff_coeff_beta_i, diff_coeff_beta_j); + + auto residual_thermal = numerics->ComputeResidual(config); + + if (ReducerStrategy) { + EdgeFluxes.SubtractBlock(iEdge, residual_thermal); + } else { + LinSysRes.SubtractBlock(iPoint, residual_thermal); + LinSysRes.AddBlock(jPoint, residual_thermal); + /* No implicit part for the preferential diffusion of heat */ + } + } +} + unsigned long CSpeciesFlameletSolver::GetEnthFromTemp(CFluidModel* fluid_model, su2double const val_temp, const su2double* scalar_solution, su2double* val_enth) { /*--- convergence criterion for temperature in [K], high accuracy needed for restarts. ---*/ diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index f5f5ed16603..05d43f13ff8 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -314,8 +314,8 @@ void CSpeciesSolver::Preprocessing(CGeometry* geometry, CSolver** solver_contain CommonPreprocessing(geometry, config, Output); } -void CSpeciesSolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { +void CSpeciesSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { /*--- Define an object to set solver specific numerics contribution. ---*/ auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) { /*--- Mass diffusivity coefficients. ---*/ diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index b67e1b1af96..af89924ea2e 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -263,8 +263,8 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai } -void CTransLMSolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { +void CTransLMSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { /*--- Define an object to set solver specific numerics contribution. ---*/ diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index f3a6d57b589..d84b56e1dc6 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -291,8 +291,8 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain AD::EndNoSharedReading(); } -void CTurbSASolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { +void CTurbSASolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { /*--- Define an object to set solver specific numerics contribution. ---*/ auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) { diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 1973b70a331..ea316e9a0cd 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -276,8 +276,8 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai AD::EndNoSharedReading(); } -void CTurbSSTSolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { +void CTurbSSTSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { /*--- Define an object to set solver specific numerics contribution. ---*/ auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) { diff --git a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp index da10ba21bd8..3dd1b222750 100644 --- a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp @@ -53,4 +53,9 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf, source_scalar.resize(nPoint, config->GetNScalars()) = su2double(0.0); lookup_scalar.resize(nPoint, config->GetNLookups()) = su2double(0.0); table_misses.resize(nPoint) = 0; + + if (config->GetPreferentialDiffusion()) { + 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); + } } diff --git a/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg b/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg new file mode 100644 index 00000000000..3f294865a17 --- /dev/null +++ b/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg @@ -0,0 +1,120 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Laminar premixed hydrogen flame with heat exchanger % +% Author: Evert Bunschoten % +% Institution: Delft University of Technology % +% Date: 01/11/2023 % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER = INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +MATH_PROBLEM= DIRECT +RESTART_SOL =YES +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= VARIABLE +INC_ENERGY_EQUATION = YES +INC_VELOCITY_INIT= (1.13, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 300.0 +THERMODYNAMIC_PRESSURE= 101325 +INC_NONDIM= DIMENSIONAL + +% -------------------- FLUID MODEL --------------------------------------- % +% +FLUID_MODEL= FLUID_FLAMELET +VISCOSITY_MODEL= FLAMELET +CONDUCTIVITY_MODEL= FLAMELET +DIFFUSIVITY_MODEL= FLAMELET +KIND_SCALAR_MODEL= FLAMELET +INTERPOLATION_METHOD= MLP +FILENAMES_INTERPOLATOR= (MLP_TD.mlp, MLP_PD.mlp, MLP_PPV.mlp, MLP_null.mlp) +PREFERENTIAL_DIFFUSION= YES + +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +% Using an artificial spark to ignite the solution at some location and iteration +FLAME_INIT_METHOD= SPARK +% Spark parameters in order: +% x-location of spark center (m) +% y-location of spark center (m) +% z-location of spark center (m) +% Spark radius (m) +% Iteration at which the artificial spark starts +% Spark iteration duration +SPARK_INIT= (0.001, 0.0004, 0, 5e-4, 10, 10) + +% Controlling variable source terms applied within the spark sphere for the spark +% duration. +SPARK_REACTION_RATES= (1000, 0, 0) + +SPECIES_INIT = (-0.49904325357252965, 2226.901776784524, 0.01446751783896619) + +% Passive reactants in flamelet problem + +CONTROLLING_VARIABLE_NAMES= (ProgressVariable, EnthalpyTot, MixtureFraction) +CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL, NULL) + +SPECIES_CLIPPING=YES +SPECIES_CLIPPING_MAX=+4.5623852432084366e-01 +8.6731375409855954e+06 1.0 +SPECIES_CLIPPING_MIN= -6.8059708053507162e-01 -4.9308262569627967e+06 0.0 + +MARKER_INLET_SPECIES = (inlet, -0.49904325357252965, 2226.901776784524, 0.01446751783896619) + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= (burner_wall, 300, hex_wall, 300) +MARKER_SYM= (sides) +INC_INLET_TYPE= VELOCITY_INLET +MARKER_INLET = (inlet, 300.000, 0.575, 1.0, 0.0, 0.0) +INC_OUTLET_TYPE= PRESSURE_OUTLET +MARKER_OUTLET= (outlet, 0.0) +MARKER_ANALYZE_AVERAGE = AREA + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER =50 +CFL_ADAPT= NO +ITER=5 +OUTPUT_WRT_FREQ= 20 +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-4 +LINEAR_SOLVER_ITER=20 +% +% -------------------- FLOW AND SPECIES NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR +MUSCL_FLOW= YES +MUSCL_SPECIES= YES +SLOPE_LIMITER_FLOW = NONE +SLOPE_LIMITER_SPECIES= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +TIME_DISCRE_SPECIES= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD = RMS_EnthalpyTot +CONV_RESIDUAL_MINVAL= -5 +CONV_STARTITER= 20 +SCREEN_OUTPUT = INNER_ITER RMS_PRESSURE RMS_ProgressVariable RMS_EnthalpyTot RMS_MixtureFraction +HISTORY_OUTPUT = WALL_TIME RMS_RES +VOLUME_OUTPUT = SOLUTION + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME = 2Dhex_BL.su2 +OUTPUT_FILES = (RESTART,PARAVIEW) +TABULAR_FORMAT = CSV +CONV_FILENAME= history +VOLUME_FILENAME= flow + + diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 51b18e48741..b891e7d9a84 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -69,6 +69,15 @@ def main(): cfd_flamelet_ch4_partial_premix.new_output = True test_list.append(cfd_flamelet_ch4_partial_premix) + # 2D planar laminar premixed hydrogen flame on isothermal burner with heat exchanger emulator (restart) + cfd_flamelet_h2 = TestCase('cfd_flamelet_h2') + cfd_flamelet_h2.cfg_dir = "flamelet/07_laminar_premixed_h2_flame_cfd" + cfd_flamelet_h2.cfg_file = "laminar_premixed_h2_flame_cfd.cfg" + cfd_flamelet_h2.test_iter = 5 + cfd_flamelet_h2.test_vals = [-10.003106, -9.843748, -3.289857, -11.338273] + cfd_flamelet_h2.new_output = True + test_list.append(cfd_flamelet_h2) + ######################### ## NEMO solver ### ######################### @@ -1045,9 +1054,9 @@ def main(): ###################################### # Aachen Turbine restart - Aachen_3D_restart = TestCase('aachen_turbine_restart') - Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" - Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" + Aachen_3D_restart = TestCase('aachen_turbine_restart') + Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" + Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 Aachen_3D_restart.test_vals = [-15.329206, -15.008622, -15.078888, -13.841072, -12.727840, -9.975729] test_list.append(Aachen_3D_restart) @@ -1065,8 +1074,8 @@ def main(): axial_stage2D.cfg_dir = "turbomachinery/axial_stage_2D" axial_stage2D.cfg_file = "Axial_stage2D.cfg" axial_stage2D.test_iter = 20 - axial_stage2D.test_vals = [0.983754, 1.534455, -2.888523, 2.606770, -2.418403, 3.087203, 106380, 106380, 5.7325, 64.711] - axial_stage2D.test_vals_aarch64 = [0.983754, 1.534455, -2.888523, 2.606770, -2.418403, 3.087203, 106380, 106380, 5.7325, 64.711] + axial_stage2D.test_vals = [0.983754, 1.534455, -2.888523, 2.606770, -2.418403, 3.087203, 106380, 106380, 5.7325, 64.711] + axial_stage2D.test_vals_aarch64 = [0.983754, 1.534455, -2.888523, 2.606770, -2.418403, 3.087203, 106380, 106380, 5.7325, 64.711] test_list.append(axial_stage2D) # 2D transonic stator restart @@ -1074,7 +1083,7 @@ def main(): transonic_stator_restart.cfg_dir = "turbomachinery/transonic_stator_2D" transonic_stator_restart.cfg_file = "transonic_stator_restart.cfg" transonic_stator_restart.test_iter = 20 - transonic_stator_restart.test_vals = [-5.011834, -3.091110, -2.757795, 1.087934, -3.544707, 2.166101, -471630, 94.868, -0.035888] + transonic_stator_restart.test_vals = [-5.011834, -3.091110, -2.757795, 1.087934, -3.544707, 2.166101, -471630, 94.868, -0.035888] transonic_stator_restart.test_vals_aarch64 = [-5.011834, -3.091110, -2.757795, 1.087934, -3.544707, 2.166101, -471630, 94.868, -0.035888] test_list.append(transonic_stator_restart) @@ -1378,7 +1387,7 @@ def main(): pywrapper_unsteadyFSI.test_iter = 4 pywrapper_unsteadyFSI.test_vals = [0, 31, 5, 58, -1.756780, -2.828276, -7.652558, -6.863929, 1.5618e-04] pywrapper_unsteadyFSI.command = TestCase.Command("mpirun -np 2", "python", "run.py") - pywrapper_unsteadyFSI.unsteady = True + pywrapper_unsteadyFSI.unsteady = True pywrapper_unsteadyFSI.multizone = True test_list.append(pywrapper_unsteadyFSI) @@ -1570,13 +1579,13 @@ def main(): ##################### # CGNS writer - cgns_writer = TestCase('cgns_writer') - cgns_writer.cfg_dir = "cgns_writer" - cgns_writer.cfg_file = "config.cfg" - cgns_writer.test_iter = 1 - cgns_writer.test_vals = [-2.974473, 0.665204, 5.068846, -7.003873] - cgns_writer.command = TestCase.Command("mpirun -n 2", "SU2_CFD") - cgns_writer.new_output = True + cgns_writer = TestCase('cgns_writer') + cgns_writer.cfg_dir = "cgns_writer" + cgns_writer.cfg_file = "config.cfg" + cgns_writer.test_iter = 1 + cgns_writer.test_vals = [-2.974473, 0.665204, 5.068846, -7.003873] + cgns_writer.command = TestCase.Command("mpirun -n 2", "SU2_CFD") + cgns_writer.new_output = True test_list.append(cgns_writer) ###################################### @@ -1806,14 +1815,14 @@ def main(): test_list.append(sphere_ffd_def_bspline) # Inviscid NACA0012 (triangles) - naca0012_cst = TestCase('naca0012_cst') - naca0012_cst.cfg_dir = "deformation/cst" - naca0012_cst.cfg_file = "naca0012.cfg" + naca0012_cst = TestCase('naca0012_cst') + naca0012_cst.cfg_dir = "deformation/cst" + naca0012_cst.cfg_file = "naca0012.cfg" naca0012_cst.test_iter = 10 naca0012_cst.test_vals = [0.000385514] #residual - naca0012_cst.command = TestCase.Command("mpirun -n 2", "SU2_DEF") - naca0012_cst.timeout = 1600 - naca0012_cst.tol = 1e-8 + naca0012_cst.command = TestCase.Command("mpirun -n 2", "SU2_DEF") + naca0012_cst.timeout = 1600 + naca0012_cst.tol = 1e-8 pass_list.append(naca0012_cst.run_def()) test_list.append(naca0012_cst) diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index 849a4f85c26..f63fb81c990 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -120,6 +120,16 @@ def main(): sudo_tutorial.command = TestCase.Command("mpirun -n 2", "SU2_CFD") test_list.append(sudo_tutorial) + ### Incompressible Combustion + + # Pre-mixed, laminar hydrogen flame with heat loss + premixed_hydrogen = TestCase('premixed_hydrogen') + premixed_hydrogen.cfg_dir = "../Tutorials/incompressible_flow/Inc_Combustion/1__premixed_hydrogen" + premixed_hydrogen.cfg_file = "H2_burner.cfg" + premixed_hydrogen.test_iter = 10 + premixed_hydrogen.test_vals = [-9.905856, -10.449512, -11.035999, -4.331440, -11.882740] + test_list.append(premixed_hydrogen) + ### Compressible Flow # Inviscid Bump diff --git a/config_template.cfg b/config_template.cfg index bf3e6bae907..74832ef98db 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -869,6 +869,11 @@ SPECIES_CLIPPING_MAX= 1.0 SPECIES_CLIPPING_MIN= 0.0 % --------------------- FLAMELET MODEL -----------------------------% +% The flamelet model uses the prompts INTERPOLATION_METHOD and FILENAMES_INTERPOLATOR +% for definition of the flamelet manifold. Either LUT or MLP can be used as options. +% If the terms "Beta_ProgVar", "Beta_Enth_Thermal", "Beta_Enth", and "Beta_Mixfrac" are +% found in the variables list of the manifold, preferential diffusion is assumed. + % % Names of the user defined (auxiliary) transport equations. USER_SCALAR_NAMES= (Y-CO) @@ -894,8 +899,27 @@ CONTROLLING_VARIABLE_NAMES= (ProgressVariable, EnthalpyTot) % controlling variables without source terms, use "zero" or "null", similar to the % option USER_SOURCE_NAMES. CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL) -% -% flamelet initialization + +% Method used to ignite the solution: +% FLAME_FRONT : the flame is initialized using a plane, defined by a point and a normal. On one side, the solution is initialized +% using 'burnt' conditions and on the other side 'unburnt' conditions. The normal points in the direction of the 'burnt' +% condition. +% SPARK : the solution is ignited through application of a set of source terms within a specified region for a set number +% of solver iterations. This artificial spark can be applied after a certain number of iterations has passed, allowing for +% the flow to evolve before igniting the mixture. +% NONE : no artificial solution ignition. +% By default, this option is set to NONE +FLAME_INIT_METHOD= FLAME_FRONT + +% Enable preferential diffusion for progress variable, total enthalpy, mixture fraction FGM problems. +% If 'YES', the preferential diffusion model from Efimov et al.(2021) is used for computing the viscous residual terms. +% If enabled, make sure the variables "Beta_ProgVar", "Beta_Enth_Thermal", "Beta_Enth", and "Beta_MixFrac", corresponding +% to the preferential diffusion terms of the progress variable, temperature, total enthalpy, and mixture fraction +% respectively are included in the manifold. +% By default, this option is disabled. +PREFERENTIAL_DIFFUSION= NO + +% FLAME_FRONT initialization % the flame is initialized using a plane, defined by a point and a normal. On one side, the solution is initialized % using 'burnt' conditions and on the other side 'unburnt' conditions. The normal points in the direction of the 'burnt' % condition. @@ -904,6 +928,21 @@ CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL) % (x7) = thickness of the reaction zone, this is the transition from unburnt to burnt conditions. % (x8) = Thickness of the 'burnt' zone, after this length, the conditions will be 'unburnt' again. FLAME_INIT= (0.004, 0.0, 0.0, 1.0, 0.0, 0.0, 0.2e-3, 1.0) + +% SPARK initialization +% the solution is ignited through application of a set of source terms within a specified region for a set number +% of solver iterations. This artificial spark can be applied after a certain number of iterations has passed, allowing for +% the flow to evolve before igniting the mixture. +% (x1,x2,x3) = spark center location. +% (x4) = spark radius where within the artificial spark is applied. +% (x5,x6) = spark iteration start and number of iterations in which the artifical spark is applied. For single-zone problems +% the number of iterations are inner-loop iterations, for multi-zone problems, outer-loop iterations are considered. +SPARK_INIT= (0.004, 0.0, 0.0, 1e-4, 100, 10) + +% Source terms (density times species time rate of change) applied to the species equations in the spark region for the duration of the spark. +% The number of terms should equate the total number of species in the flamelet problem. +SPARK_REACTION_RATES= (1000, 0, 0) + % % --------------------- INVERSE DESIGN SIMULATION -----------------------------% % diff --git a/meson_scripts/init.py b/meson_scripts/init.py index 787bfb56e2f..bc58f9a9564 100755 --- a/meson_scripts/init.py +++ b/meson_scripts/init.py @@ -71,7 +71,7 @@ def init_submodules( github_repo_coolprop = "https://github.com/CoolProp/CoolProp" sha_version_mel = "46205ab019e5224559091375a6d71aabae6bc5b9" github_repo_mel = "https://github.com/pcarruscag/MEL" - sha_version_mlpcpp = "665c45b7d3533c977eb1f637918d5b8b75c07d3b" + sha_version_mlpcpp = "c19c53ea2b85ccfb185f1c6c87044dc0b5bc7ae0" github_repo_mlpcpp = "https://github.com/EvertBunschoten/MLPCpp" medi_name = "MeDiPack" diff --git a/subprojects/MLPCpp b/subprojects/MLPCpp index 665c45b7d35..c19c53ea2b8 160000 --- a/subprojects/MLPCpp +++ b/subprojects/MLPCpp @@ -1 +1 @@ -Subproject commit 665c45b7d3533c977eb1f637918d5b8b75c07d3b +Subproject commit c19c53ea2b85ccfb185f1c6c87044dc0b5bc7ae0