diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index f3b1abd98764..ccaba0ac2471 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -379,6 +379,16 @@ class CConfig { su2double *Surface_PressureDrop; /*!< \brief Pressure drop between boundaries. */ su2double* Surface_Species_0; /*!< \brief Average Species_0 at the boundaries. */ su2double* Surface_Species_Variance; /*!< \brief Species Variance at the boundaries. */ + su2double* Surface_Scalar_00; /*!< \brief Average of scalar 0 at the boundaries. */ + su2double* Surface_Scalar_01; /*!< \brief Average of scalar 1 at the boundaries. */ + su2double* Surface_Scalar_02; /*!< \brief Average of scalar 2 at the boundaries. */ + su2double* Surface_Scalar_03; /*!< \brief Average of scalar 3 at the boundaries. */ + su2double* Surface_Scalar_04; /*!< \brief Average of scalar 4 at the boundaries. */ + su2double* Surface_Scalar_05; /*!< \brief Average of scalar 5 at the boundaries. */ + su2double* Surface_Scalar_06; /*!< \brief Average of scalar 6 at the boundaries. */ + su2double* Surface_Scalar_07; /*!< \brief Average of scalar 7 at the boundaries. */ + su2double* Surface_Scalar_08; /*!< \brief Average of scalar 8 at the boundaries. */ + su2double* Surface_Scalar_09; /*!< \brief Average of scalar 9 at the boundaries. */ su2double *Surface_DC60; /*!< \brief Specified surface DC60 for nacelle boundaries. */ su2double *Surface_IDC; /*!< \brief Specified IDC for nacelle boundaries. */ su2double *Surface_IDC_Mach; /*!< \brief Specified IDC mach for nacelle boundaries. */ @@ -714,6 +724,7 @@ class CConfig { *Marker_WallFunctions, /*!< \brief Markers for which wall functions must be applied. */ *Marker_SobolevBC; /*!< \brief Markers in the gradient solver */ + bool initial_PyCustom; /*!< \brief flag for using custom python boundary conditions */ unsigned short nConfig_Files; /*!< \brief Number of config files for multiphysics problems. */ string *Config_Filenames; /*!< \brief List of names for configuration files. */ SST_OPTIONS *SST_Options; /*!< \brief List of modifications/corrections/versions of SST turbulence model.*/ @@ -767,6 +778,7 @@ class CConfig { unsigned short ActDisk_Jump; /*!< \brief Format of the output files. */ unsigned long StartWindowIteration; /*!< \brief Starting Iteration for long time Windowing apporach . */ unsigned short nCFL_AdaptParam; /*!< \brief Number of CFL parameters provided in config. */ + bool Initial_All_PyCustom; /*!< \brief Python customizable initial condition */ bool CFL_Adapt; /*!< \brief Use adaptive CFL number. */ bool HB_Precondition; /*!< \brief Flag to turn on harmonic balance source term preconditioning */ su2double RefArea, /*!< \brief Reference area for coefficient computation. */ @@ -1215,6 +1227,27 @@ class CConfig { su2double* Species_Init; /*!< \brief Initial uniform value for scalar transport. */ unsigned short nSpecies_Init; /*!< \brief Number of entries of SPECIES_INIT */ + /*--- flamelet subsolver ---*/ + su2double flame_thickness; + su2double flame_burnt_thickness; + su2double flame_offset[3]; + su2double flame_normal[3]; + + /*--- lookup table ---*/ + unsigned short n_scalars; /* number of transported scalars for the flamelet LUT approach*/ + unsigned short n_lookups; /* number of lookud up variables */ + unsigned short n_table_sources; /* the number of transported scalar source terms for the LUT */ + unsigned short n_user_scalars; + unsigned short n_user_sources; + unsigned short n_control_vars; + + vector table_scalar_names; /*!< \brief vector to store names of scalar variables. */ + vector table_source_names; /*!< \brief vector to store names of scalar source variables. */ + string* table_lookup_names; /*!< \brief vector to store names of look up variables. */ + string file_name_lut; /*!< \brief file name of the look up table. */ + string* user_scalar_names; + string* user_source_names; + /*! * \brief Set the default values of config options not set in the config file using another config object. * \param config - Config object to use the default values from. @@ -2106,6 +2139,110 @@ class CConfig { */ bool GetSpecies_StrongBC() const { return Species_StrongBC; } + /*! + * \brief Get the flame offset for flamelet model initialization + * \return flame offset for flamelet model initialization + */ + su2double *GetFlameOffset(void) { return flame_offset; } + + /*! + * \brief Get the flame normal for flamelet model initialization + * \return flame offset for flamelet model initialization + */ + su2double *GetFlameNormal(void) { return flame_normal; } + + /*! + * \brief Get the flame thickness for flamelet model initialization + * \return flame thickness for flamelet model initialization + */ + su2double GetFlameThickness(void) { return flame_thickness; } + + /*! + * \brief Get the burnt region thickness for flamelet mdoel initialization + * \return flame thickness for flamelet model initialization + */ + su2double GetFlameBurntThickness(void) { return flame_burnt_thickness; } + + /*! + * \brief Set the number of scalars for flamelet model. + */ + void SetNScalars(unsigned short n_scalars) { this->n_scalars = n_scalars; } + + /*! + * \brief Set the number of controlling variables for flamelet model. + */ + void SetNControlVars(unsigned short n_control_vars) { this->n_control_vars = n_control_vars; } + + /*! + * \brief Get the number of control variables for flamelet model. + */ + unsigned short GetNControlVars(void) const { return n_control_vars; } + + /*! + * \brief Get the number of transported scalars for flamelet model. + */ + unsigned short GetNScalars(void) const { return n_scalars; } + + /*! + * \brief Get the number of user scalars for flamelet model. + */ + unsigned short GetNUserScalars(void) const { return n_user_scalars; } + + /*! + * \brief Get the name of the user scalar. + */ + string GetUserScalarName(unsigned short i_user_scalar) const { if(n_user_scalars > 0) return user_scalar_names[i_user_scalar]; else return "NONE"; } + + /*! + * \brief Get the name of the user scalar source term. + */ + string GetUserSourceName(unsigned short i_user_source) const { if(n_user_sources > 0) return user_source_names[i_user_source]; else return "NONE"; } + + /*! + * \brief Get the number of transported scalars for combustion + */ + unsigned short GetNLookups(void) const { return n_lookups; } + + void SetNLUTSources(unsigned short n_table_sources) { this->n_table_sources = n_table_sources; } + + /*! + * \brief Get the number of transported scalars source terms for combustion + */ + unsigned short GetNLUTSources(void) const { return n_table_sources; } + + /*! + * \brief Store the names of scalar variables that are being solved + * \param[out] stores the names in vector table_scalar_names + */ + inline void SetLUTScalarNames(vector &table_scalar_names) { this->table_scalar_names = table_scalar_names; } + + /*! + * \brief Get the name of the independent variable from the lookup table + */ + string GetLUTScalarName(unsigned short i_scalar) const { return table_scalar_names[i_scalar]; } + + /*! + * \brief Get the name of the variable that we want to retrieve from the lookup table + */ + string GetLUTLookupName(unsigned short i_lookup) const { return table_lookup_names[i_lookup]; } + + /*! + * \brief Store the names of scalar source term variables + * \param[out] stores the names in vector table_source_names + */ + inline void SetLUTSourceNames(vector &table_source_names) { this->table_source_names = table_source_names; } + + /*! + * \brief Get the scalar source term name i_source + */ + string GetLUTSourceName(unsigned short i_source) const { return table_source_names[i_source]; } + + /*! + * \brief Get the file name of the look up table + * \return File name of the look up table + */ + string GetFileNameLUT(void){ return file_name_lut; }; + /*! * \brief Get the Young's modulus of elasticity. * \return Value of the Young's modulus of elasticity. @@ -2113,15 +2250,15 @@ class CConfig { su2double GetElasticyMod(unsigned short id_val) const { return ElasticityMod[id_val]; } /*! - * \brief Decide whether to apply DE effects to the model. - * \return TRUE if the DE effects are to be applied, FALSE otherwise. - */ + * \brief Decide whether to apply DE effects to the model. + * \return TRUE if the DE effects are to be applied, FALSE otherwise. + */ bool GetDE_Effects(void) const { return DE_Effects; } /*! - * \brief Decide whether to predict the DE effects for the next time step. - * \return TRUE if the DE effects are to be applied, FALSE otherwise. - */ + * \brief Decide whether to predict the DE effects for the next time step. + * \return TRUE if the DE effects are to be applied, FALSE otherwise. + */ bool GetDE_Predicted(void); /*! @@ -3022,6 +3159,12 @@ class CConfig { */ unsigned short GetnMarker_PyCustom(void) const { return nMarker_PyCustom; } + /*! + * \brief Get Python customizable initial condition. + * \return true if customizable initial condition + */ + bool GetInitial_PyCustom(void) const { return initial_PyCustom; } + /*! * \brief Get the total number of moving markers. * \return Total number of moving markers. @@ -3414,6 +3557,13 @@ class CConfig { */ void SetMarker_All_PyCustom(unsigned short val_marker, unsigned short val_PyCustom) { Marker_All_PyCustom[val_marker] = val_PyCustom; } + /*! + * \brief Set if the initial condition is going to be customized in Python val_PyCustom + * (read from the config file). + * \param[in] val_PyCustom - 0 or 1 depending if the the initial condition is going to be customized in Python. + */ + void SetInitial_All_PyCustom(unsigned short val_PyCustom) { Initial_All_PyCustom = val_PyCustom; } + /*! * \brief Set if a marker val_marker is going to be periodic val_perbound * (read from the config file). @@ -3564,6 +3714,12 @@ class CConfig { */ unsigned short GetMarker_All_PyCustom(unsigned short val_marker) const { return Marker_All_PyCustom[val_marker];} + /*! + * \brief Get the Python customization for the initial condition + * \return 0 or 1 depending if the initial condition is going to be customized in Python. + */ + unsigned short GetInitial_All_PyCustom(unsigned short val_initial) const { return Initial_All_PyCustom;} + /*! * \brief Get the airfoil sections in the slicing process. * \param[in] val_section - Index of the section. @@ -7727,6 +7883,22 @@ class CConfig { */ void SetSurface_Species_Variance(unsigned short val_marker, su2double val_surface_species_variance) { Surface_Species_Variance[val_marker] = val_surface_species_variance; } + /*! + * \brief Set the average of scalar_0 at the surface. + * \param[in] val_marker - Index corresponding to boundary. + * \param[in] val_surface_scalar_0 - Value of avg species_0. + */ + void SetSurface_Scalar_00(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_00[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_01(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_01[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_02(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_02[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_03(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_03[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_04(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_04[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_05(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_05[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_06(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_06[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_07(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_07[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_08(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_08[val_marker] = val_surface_scalar; } + void SetSurface_Scalar_09(unsigned short val_marker, su2double val_surface_scalar) { Surface_Scalar_09[val_marker] = val_surface_scalar; } + /*! * \brief Get the back pressure (static) at an outlet boundary. * \param[in] val_index - Index corresponding to the outlet boundary. @@ -8007,6 +8179,22 @@ class CConfig { */ su2double GetSurface_Species_Variance(unsigned short val_marker) const { return Surface_Species_Variance[val_marker]; } + /*! + * \brief Get avg scalar_0 at a boundary. + * \param[in] val_index - Index corresponding to the boundary. + * \return The avg species_0. + */ + su2double GetSurface_Scalar_00(unsigned short val_marker) const { return Surface_Scalar_00[val_marker]; } + su2double GetSurface_Scalar_01(unsigned short val_marker) const { return Surface_Scalar_01[val_marker]; } + su2double GetSurface_Scalar_02(unsigned short val_marker) const { return Surface_Scalar_02[val_marker]; } + su2double GetSurface_Scalar_03(unsigned short val_marker) const { return Surface_Scalar_03[val_marker]; } + su2double GetSurface_Scalar_04(unsigned short val_marker) const { return Surface_Scalar_04[val_marker]; } + su2double GetSurface_Scalar_05(unsigned short val_marker) const { return Surface_Scalar_05[val_marker]; } + su2double GetSurface_Scalar_06(unsigned short val_marker) const { return Surface_Scalar_06[val_marker]; } + su2double GetSurface_Scalar_07(unsigned short val_marker) const { return Surface_Scalar_07[val_marker]; } + su2double GetSurface_Scalar_08(unsigned short val_marker) const { return Surface_Scalar_08[val_marker]; } + su2double GetSurface_Scalar_09(unsigned short val_marker) const { return Surface_Scalar_09[val_marker]; } + /*! * \brief Get the back pressure (static) at an outlet boundary. * \param[in] val_index - Index corresponding to the outlet boundary. diff --git a/Common/include/containers/CTrapezoidalMap.hpp b/Common/include/containers/CTrapezoidalMap.hpp index 982f94309926..afbeb2d86ea1 100644 --- a/Common/include/containers/CTrapezoidalMap.hpp +++ b/Common/include/containers/CTrapezoidalMap.hpp @@ -29,9 +29,10 @@ #include #include +#include -#include "../../Common/include/linear_algebra/blas_structure.hpp" #include "../../Common/include/toolboxes/CSquareMatrixCM.hpp" +#include "../../Common/include/linear_algebra/blas_structure.hpp" /*! * \class CTrapezoidalMap diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index c4e5067245c2..0fc0128ade61 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -1275,6 +1275,13 @@ class CGeometry { */ void UpdateCustomBoundaryConditions(CGeometry **geometry_container, CConfig *config); + /*! + * \brief Update the multi-grid structure for the customized initial conditions + * \param geometry_container - Geometrical definition. + * \param config - Definition of the particular problem. + */ + void UpdateCustomInitialConditions(CGeometry **geometry_container, CConfig *config); + /*! * \brief A virtual member. * \param config - Config diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 72347466647c..a6d47f5e8a77 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -471,28 +471,44 @@ enum RUNTIME_TYPE { RUNTIME_ADJSPECIES_SYS = 26,/*!< \brief One-physics case, the code is solving the adjoint species model. */ }; -const int FLOW_SOL = 0; /*!< \brief Position of the mean flow solution in the solver container array. */ -const int ADJFLOW_SOL = 1; /*!< \brief Position of the continuous adjoint flow solution in the solver container array. */ - -const int TURB_SOL = 2; /*!< \brief Position of the turbulence model solution in the solver container array. */ -const int ADJTURB_SOL = 3; /*!< \brief Position of the continuous adjoint turbulence solution in the solver container array. */ - -const int TRANS_SOL = 4; /*!< \brief Position of the transition model solution in the solver container array. */ -const int HEAT_SOL = 5; /*!< \brief Position of the heat equation in the solution solver array. */ -const int ADJHEAT_SOL = 6; /*!< \brief Position of the adjoint heat equation in the solution solver array. */ -const int RAD_SOL = 7; /*!< \brief Position of the radiation equation in the solution solver array. */ -const int ADJRAD_SOL = 8; /*!< \brief Position of the continuous adjoint turbulence solution in the solver container array. */ - -const int MESH_SOL = 9; /*!< \brief Position of the mesh solver. */ -const int ADJMESH_SOL = 10; /*!< \brief Position of the adjoint of the mesh solver. */ - -const int SPECIES_SOL = 11; /*!< \brief Position of the species solver. */ -const int ADJSPECIES_SOL = 12; /*!< \brief Position of the adjoint of the species solver. */ - -const int FEA_SOL = 0; /*!< \brief Position of the FEA equation in the solution solver array. */ -const int ADJFEA_SOL = 1; /*!< \brief Position of the FEA adjoint equation in the solution solver array. */ - -const int TEMPLATE_SOL = 0; /*!< \brief Position of the template solution. */ + // nijso: the values of int are important for some cases (fsi) + enum SOLVER_TYPE : const int { + FLOW_SOL=0, /*!< \brief Position of the mean flow solution in the solver container array. */ + ADJFLOW_SOL=1, /*!< \brief Position of the continuous adjoint flow solution in the solver container array. */ + TURB_SOL=2, /*!< \brief Position of the turbulence model solution in the solver container array. */ + ADJTURB_SOL=3, /*!< \brief Position of the continuous adjoint turbulence solution in the solver container array. */ + TRANS_SOL=4, /*!< \brief Position of the transition model solution in the solver container array. */ + HEAT_SOL=5, /*!< \brief Position of the heat equation in the solution solver array. */ + ADJHEAT_SOL=6, /*!< \brief Position of the adjoint heat equation in the solution solver array. */ + RAD_SOL=7, /*!< \brief Position of the radiation equation in the solution solver array. */ + ADJRAD_SOL=8, /*!< \brief Position of the continuous adjoint turbulence solution in the solver container array. */ + MESH_SOL=9, /*!< \brief Position of the mesh solver. */ + ADJMESH_SOL=10, /*!< \brief Position of the adjoint of the mesh solver. */ + SPECIES_SOL=11, /*!< \brief Position of the species solver. */ + ADJSPECIES_SOL=12, /*!< \brief Position of the adjoint of the species solver. */ + FEA_SOL=0, /*!< \brief Position of the Finite Element flow solution in the solver container array. */ + ADJFEA_SOL=1, /*!< \brief Position of the continuous adjoint Finite Element flow solution in the solver container array. */ + TEMPLATE_SOL=0, /*!< \brief Position of the template solution. */ + }; + + static const MapType SolverType_Map = { + MakePair("FLOW_SOL", FLOW_SOL) + MakePair("ADJFLOW_SOL", ADJFLOW_SOL) + MakePair("TURB_SOL", TURB_SOL) + MakePair("ADJTURB_SOL", ADJTURB_SOL) + MakePair("FEA_SOL", FEA_SOL) + MakePair("ADJFEA_SOL", ADJFEA_SOL) + MakePair("TEMPLATE_SOL", TEMPLATE_SOL) + MakePair("TRANS_SOL", TRANS_SOL) + MakePair("HEAT_SOL", HEAT_SOL) + MakePair("ADJHEAT_SOL", ADJHEAT_SOL) + MakePair("RAD_SOL", RAD_SOL) + MakePair("ADJRAD_SOL", ADJRAD_SOL) + MakePair("MESH_SOL", MESH_SOL) + MakePair("ADJMESH_SOL", ADJMESH_SOL) + MakePair("SPECIES_SOL", SPECIES_SOL) + MakePair("ADJSPECIES_SOL", ADJSPECIES_SOL) + }; const int CONV_TERM = 0; /*!< \brief Position of the convective terms in the numerics container array. */ const int VISC_TERM = 1; /*!< \brief Position of the viscous terms in the numerics container array. */ @@ -550,6 +566,7 @@ enum ENUM_FLUIDMODEL { SU2_NONEQ = 8, /*!< \brief User defined gas model for nonequilibrium flow. */ FLUID_MIXTURE = 9, /*!< \brief Species mixture model. */ COOLPROP = 10, /*!< \brief Thermodynamics library. */ + FLUID_FLAMELET = 11, /*!< \brief lookup table (LUT) method for premixed flamelets. */ }; static const MapType FluidModel_Map = { MakePair("STANDARD_AIR", STANDARD_AIR) @@ -563,6 +580,7 @@ static const MapType FluidModel_Map = { MakePair("SU2_NONEQ", SU2_NONEQ) MakePair("FLUID_MIXTURE", FLUID_MIXTURE) MakePair("COOLPROP", COOLPROP) + MakePair("FLUID_FLAMELET", FLUID_FLAMELET) }; /*! @@ -652,12 +670,14 @@ enum class VISCOSITYMODEL { CONSTANT, /*!< \brief Constant viscosity. */ SUTHERLAND, /*!< \brief Sutherlands Law viscosity. */ POLYNOMIAL, /*!< \brief Polynomial viscosity. */ + FLAMELET, /*!< \brief LUT method for flamelets */ COOLPROP, /*!< \brief CoolProp viscosity. */ }; static const MapType ViscosityModel_Map = { MakePair("CONSTANT_VISCOSITY", VISCOSITYMODEL::CONSTANT) MakePair("SUTHERLAND", VISCOSITYMODEL::SUTHERLAND) MakePair("POLYNOMIAL_VISCOSITY", VISCOSITYMODEL::POLYNOMIAL) + MakePair("FLAMELET", VISCOSITYMODEL::FLAMELET) MakePair("COOLPROP", VISCOSITYMODEL::COOLPROP) }; @@ -680,12 +700,14 @@ enum class CONDUCTIVITYMODEL { CONSTANT, /*!< \brief Constant thermal conductivity. */ CONSTANT_PRANDTL, /*!< \brief Constant Prandtl number. */ POLYNOMIAL, /*!< \brief Polynomial thermal conductivity. */ + FLAMELET, /*!< \brief LUT method for flamelets */ COOLPROP, /*!< \brief COOLPROP thermal conductivity. */ }; static const MapType ConductivityModel_Map = { MakePair("CONSTANT_CONDUCTIVITY", CONDUCTIVITYMODEL::CONSTANT) MakePair("CONSTANT_PRANDTL", CONDUCTIVITYMODEL::CONSTANT_PRANDTL) MakePair("POLYNOMIAL_CONDUCTIVITY", CONDUCTIVITYMODEL::POLYNOMIAL) + MakePair("FLAMELET", CONDUCTIVITYMODEL::FLAMELET) MakePair("COOLPROP", CONDUCTIVITYMODEL::COOLPROP) }; @@ -709,6 +731,7 @@ enum class DIFFUSIVITYMODEL { CONSTANT_SCHMIDT, /*!< \brief Constant Schmidt number for mass diffusion in scalar transport. */ UNITY_LEWIS, /*!< \brief Unity Lewis model for mass diffusion in scalar transport. */ CONSTANT_LEWIS, /*!< \brief Different Lewis number model for mass diffusion in scalar transport. */ + FLAMELET, /*!< \brief flamelet model for tabulated chemistry, diffusivity from lookup table */ }; static const MapType Diffusivity_Model_Map = { @@ -716,6 +739,7 @@ static const MapType Diffusivity_Model_Map = { MakePair("CONSTANT_SCHMIDT", DIFFUSIVITYMODEL::CONSTANT_SCHMIDT) MakePair("UNITY_LEWIS", DIFFUSIVITYMODEL::UNITY_LEWIS) MakePair("CONSTANT_LEWIS", DIFFUSIVITYMODEL::CONSTANT_LEWIS) + MakePair("FLAMELET", DIFFUSIVITYMODEL::FLAMELET) }; /*! @@ -1280,11 +1304,34 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh */ enum class SPECIES_MODEL { NONE, /*!< \brief No scalar transport model. */ - SPECIES_TRANSPORT, /*!< \brief Passive scalar transport model. */ + SPECIES_TRANSPORT, /*!< \brief species transport model. */ + FLAMELET, /*!< \brief flamelet model. */ }; static const MapType Species_Model_Map = { MakePair("NONE", SPECIES_MODEL::NONE) MakePair("SPECIES_TRANSPORT", SPECIES_MODEL::SPECIES_TRANSPORT) + MakePair("FLAMELET", SPECIES_MODEL::FLAMELET) +}; + +/* the order matters: */ +/*! +* \brief progress variable and enthalpy are the first and second entry in the lookup table. +* CO and NOX are auxiliary transport equations that we solve +*/ +enum FLAMELET_SCALAR_VARIABLES { + I_PROGVAR, + I_ENTH, +}; + +/*! + * \brief the source terms for the flamelet method + */ +enum FLAMELET_SCALAR_SOURCES { + I_SRC_TOT_PROGVAR, + I_SRC_PROD_CO, + I_SRC_CONS_CO, + I_SRC_PROD_NOX, + I_SRC_CONS_NOX, }; /*! @@ -1869,6 +1916,16 @@ enum ENUM_OBJECTIVE { SURFACE_PRESSURE_DROP = 56, /*!< \brief Pressure drop objective function definition. */ SURFACE_SPECIES_0 = 58, /*!< \brief Surface Avg. Species_0 objective function definition. */ SURFACE_SPECIES_VARIANCE = 59,/*!< \brief Species Variance objective function definition. */ + SURFACE_SCALAR_00 = 80, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_01 = 81, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_02 = 82, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_03 = 83, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_04 = 84, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_05 = 85, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_06 = 86, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_07 = 87, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_08 = 88, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ + SURFACE_SCALAR_09 = 89, /*!< \brief Surface Avg. Scalar_0 objective function definition. */ CUSTOM_OBJFUNC = 31, /*!< \brief Custom objective function definition. */ FLOW_ANGLE_OUT = 46, MASS_FLOW_IN = 47, @@ -1914,6 +1971,16 @@ static const MapType Objective_Map = { MakePair("SURFACE_PRESSURE_DROP", SURFACE_PRESSURE_DROP) MakePair("SURFACE_SPECIES_0", SURFACE_SPECIES_0) MakePair("SURFACE_SPECIES_VARIANCE", SURFACE_SPECIES_VARIANCE) + MakePair("SURFACE_SCALAR_00", SURFACE_SCALAR_00) + MakePair("SURFACE_SCALAR_01", SURFACE_SCALAR_01) + MakePair("SURFACE_SCALAR_02", SURFACE_SCALAR_02) + MakePair("SURFACE_SCALAR_03", SURFACE_SCALAR_03) + MakePair("SURFACE_SCALAR_04", SURFACE_SCALAR_04) + MakePair("SURFACE_SCALAR_05", SURFACE_SCALAR_05) + MakePair("SURFACE_SCALAR_06", SURFACE_SCALAR_06) + MakePair("SURFACE_SCALAR_07", SURFACE_SCALAR_07) + MakePair("SURFACE_SCALAR_08", SURFACE_SCALAR_08) + MakePair("SURFACE_SCALAR_09", SURFACE_SCALAR_09) MakePair("CUSTOM_OBJFUNC", CUSTOM_OBJFUNC) MakePair("FLOW_ANGLE_OUT", FLOW_ANGLE_OUT) MakePair("MASS_FLOW_IN", MASS_FLOW_IN) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index f0c0d8ea91a4..2ae93dd4a331 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -830,7 +830,7 @@ void CConfig::SetPointersNull(void) { Marker_CfgFile_TurbomachineryFlag = nullptr; Marker_All_TurbomachineryFlag = nullptr; Marker_CfgFile_MixingPlaneInterface = nullptr; Marker_All_MixingPlaneInterface = nullptr; - Marker_CfgFile_PyCustom = nullptr; Marker_All_PyCustom = nullptr; + Marker_CfgFile_PyCustom = nullptr; Marker_All_PyCustom = nullptr; Initial_All_PyCustom = false; Marker_DV = nullptr; Marker_Moving = nullptr; Marker_Monitoring = nullptr; Marker_Designing = nullptr; Marker_GeoEval = nullptr; Marker_Plotting = nullptr; @@ -929,6 +929,17 @@ void CConfig::SetPointersNull(void) { Surface_DC60 = nullptr; Surface_IDC = nullptr; Surface_Species_Variance = nullptr; Surface_Species_0 = nullptr; + Surface_Scalar_00 = nullptr; + Surface_Scalar_01 = nullptr; + Surface_Scalar_02 = nullptr; + Surface_Scalar_03 = nullptr; + Surface_Scalar_04 = nullptr; + Surface_Scalar_05 = nullptr; + Surface_Scalar_06 = nullptr; + Surface_Scalar_07 = nullptr; + Surface_Scalar_08 = nullptr; + Surface_Scalar_09 = nullptr; + Outlet_MassFlow = nullptr; Outlet_Density = nullptr; Outlet_Area = nullptr; Surface_Uniformity = nullptr; Surface_SecondaryStrength = nullptr; Surface_SecondOverUniform = nullptr; @@ -1353,6 +1364,24 @@ void CConfig::SetConfig_Options() { /*!\brief SPECIES_CLIPPING \n DESCRIPTION: Use strong inlet and outlet BC in the species solver \n DEFAULT: false \ingroup Config*/ addBoolOption("SPECIES_USE_STRONG_BC", Species_StrongBC, false); + /*!\brief FLAME_OFFSET \n DESCRIPTION: Offset for flame initialization using the flamelet model \ingroup Config*/ + flame_offset[0] = 0.0; + flame_offset[1] = 0.0; + flame_offset[2] = 0.0; + addDoubleArrayOption("FLAME_OFFSET", 3,flame_offset); + + /*!\brief FLAME_THICKNESS \n DESCRIPTION: Thickness for flame initialization using the flamelet model \ingroup Config*/ + addDoubleOption("FLAME_THICKNESS", flame_thickness, 0.5e-3); + + /*!\brief FLAME_NORMAL \n DESCRIPTION: Normal for flame initialization using the flamelet model \ingroup Config*/ + flame_normal[0] = 1.0; + flame_normal[1] = 0.0; + flame_normal[2] = 0.0; + addDoubleArrayOption("FLAME_NORMAL", 3, flame_normal); + + /*!\brief FLAME_BURNT_THICKNESS \n DESCRIPTION: burnt thickness for flame initialization using the flamelet model \ingroup Config*/ + addDoubleOption("FLAME_BURNT_THICKNESS", flame_burnt_thickness, 1); + /*--- Options related to mass diffusivity and thereby the species solver. ---*/ /*!\brief DIFFUSIVITY_MODEL\n DESCRIPTION: mass diffusivity model \n DEFAULT constant disffusivity \ingroup Config*/ @@ -1489,6 +1518,10 @@ void CConfig::SetConfig_Options() { /*!\brief MARKER_PYTHON_CUSTOM\n DESCRIPTION: Python customizable marker(s) \ingroup Config*/ addStringListOption("MARKER_PYTHON_CUSTOM", nMarker_PyCustom, Marker_PyCustom); + /*!\brief INITIAL_PYTHON_CUSTOM\n DESCRIPTION: flag for using Python customizable initial condition + * \ingroup Config + */ + addBoolOption("INITIAL_PYTHON_CUSTOM", initial_PyCustom, false); /*!\brief MARKER_WALL_FUNCTIONS\n DESCRIPTION: Viscous wall markers for which wall functions must be applied. Format: (Wall function marker, wall function type, ...) \ingroup Config*/ @@ -2081,6 +2114,18 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Determine if we need to allocate memory to store the multizone residual. \n DEFAULT: true (temporarily) */ addBoolOption("MULTIZONE_RESIDUAL", Multizone_Residual, false); + /*!\brief File name of the flamelet look up table.*/ + addStringOption("FILENAME_LUT", file_name_lut, string("LUT")); + + /* DESCRIPTION: Names of the passive lookup variables for flamelet LUT */ + addStringListOption("LOOKUP_NAMES", n_lookups, table_lookup_names); + + /* DESCRIPTION: Names of the user transport equations solved in the flamelet problem. */ + addStringListOption("USER_SCALAR_NAMES", n_user_scalars, user_scalar_names); + + /* DESCRIPTION: Names of the user scalar source terms. */ + addStringListOption("USER_SOURCE_NAMES", n_user_sources, user_source_names); + /*!\brief CONV_FILENAME \n DESCRIPTION: Output file convergence history (w/o extension) \n DEFAULT: history \ingroup Config*/ addStringOption("CONV_FILENAME", Conv_FileName, string("history")); /*!\brief BREAKDOWN_FILENAME \n DESCRIPTION: Output file forces breakdown \ingroup Config*/ @@ -3344,6 +3389,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i (Kind_FluidModel == IDEAL_GAS) || (Kind_FluidModel == INC_IDEAL_GAS) || (Kind_FluidModel == FLUID_MIXTURE) || + (Kind_FluidModel == FLUID_FLAMELET) || (Kind_FluidModel == INC_IDEAL_GAS_POLY) || (Kind_FluidModel == CONSTANT_DENSITY)); bool noneq_gas = ((Kind_FluidModel == MUTATIONPP) || @@ -3659,6 +3705,16 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i case SURFACE_PRESSURE_DROP: case SURFACE_SPECIES_0: case SURFACE_SPECIES_VARIANCE: + case SURFACE_SCALAR_00: + case SURFACE_SCALAR_01: + case SURFACE_SCALAR_02: + case SURFACE_SCALAR_03: + case SURFACE_SCALAR_04: + case SURFACE_SCALAR_05: + case SURFACE_SCALAR_06: + case SURFACE_SCALAR_07: + case SURFACE_SCALAR_08: + case SURFACE_SCALAR_09: case CUSTOM_OBJFUNC: if (Kind_ObjFunc[iObj] != Obj_0) { SU2_MPI::Error("The following objectives can only be used for the first surface in a multi-objective \n" @@ -3666,7 +3722,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i "INVERSE_DESIGN_PRESSURE, INVERSE_DESIGN_HEATFLUX, THRUST_COEFFICIENT, TORQUE_COEFFICIENT\n" "FIGURE_OF_MERIT, SURFACE_TOTAL_PRESSURE, SURFACE_STATIC_PRESSURE, SURFACE_MASSFLOW\n" "SURFACE_UNIFORMITY, SURFACE_SECONDARY, SURFACE_MOM_DISTORTION, SURFACE_SECOND_OVER_UNIFORM\n" - "SURFACE_PRESSURE_DROP, SURFACE_STATIC_TEMPERATURE, SURFACE_SPECIES_0\n" + "SURFACE_PRESSURE_DROP, SURFACE_STATIC_TEMPERATURE, SURFACE_SPECIES_0, SURFACE_SCALAR_\n" "SURFACE_SPECIES_VARIANCE, CUSTOM_OBJFUNC.\n", CURRENT_FUNCTION); } break; @@ -3676,9 +3732,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } } - if (Kind_ObjFunc[0] == CUSTOM_OBJFUNC && CustomObjFunc.empty() && !Multizone_Problem) { - SU2_MPI::Error("The expression for the custom objective function was not set.\n" - "For example, CUSTOM_OBJFUNC= LIFT/DRAG", CURRENT_FUNCTION); + if (nObj > 0){ + if (Kind_ObjFunc[0] == CUSTOM_OBJFUNC && CustomObjFunc.empty() && !Multizone_Problem) { + SU2_MPI::Error("The expression for the custom objective function was not set.\n" + "For example, CUSTOM_OBJFUNC= LIFT/DRAG", CURRENT_FUNCTION); + } } /*--- Check for unsteady problem ---*/ @@ -3878,7 +3936,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } break; default: - if (nSpecies_Init + 1 != 1) SU2_MPI::Error("Viscosity model not available.", CURRENT_FUNCTION); + if (nSpecies_Init + 1 != 1) SU2_MPI::Error("Fluid mixture: viscosity model not available.", CURRENT_FUNCTION); break; } @@ -3925,6 +3983,34 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } } + + if (Kind_Species_Model == SPECIES_MODEL::FLAMELET) { + + if (Kind_FluidModel != FLUID_FLAMELET) { + SU2_MPI::Error("The use of SCALAR_MODEL= FLAMELET requires the FLUID_MODEL option to be FLUID_FLAMELET", + CURRENT_FUNCTION); + } + + if (Kind_DensityModel != INC_DENSITYMODEL::VARIABLE) { + SU2_MPI::Error("The use of FLUID_FLAMELET requires the INC_DENSITY_MODEL option to be VARIABLE", + CURRENT_FUNCTION); + } + + if (Kind_ConductivityModel != CONDUCTIVITYMODEL::FLAMELET) { + SU2_MPI::Error("The use of FLUID_FLAMELET requires the CONDUCTIVITY_MODEL option to be FLAMELET", + CURRENT_FUNCTION); + } + + if (Kind_Diffusivity_Model != DIFFUSIVITYMODEL::FLAMELET) { + SU2_MPI::Error("The use of FLUID_FLAMELET requires the DIFFUSIVITY_MODEL option to be FLAMELET", + CURRENT_FUNCTION); + } + if (Kind_ViscosityModel != VISCOSITYMODEL::FLAMELET) { + SU2_MPI::Error("The use of FLUID_FLAMELET requires the VISCOSITY_MODEL option to be FLAMELET", + CURRENT_FUNCTION); + } + } + /*--- Check for Measurement System ---*/ if (SystemMeasurements == US && !standard_air) { @@ -4877,9 +4963,9 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i Kind_FluidModel = CONSTANT_DENSITY; } - /*--- Energy equation must be active for any fluid models other than constant density. ---*/ - if ((Kind_DensityModel != INC_DENSITYMODEL::CONSTANT) && (Kind_Species_Model==SPECIES_MODEL::NONE)) Energy_Equation = true; + /*--- For the flamelet combustion model, energy equation is a passive field, we lookup T and write it to the field ---*/ + if (Kind_Species_Model == SPECIES_MODEL::FLAMELET ) Energy_Equation = false; if (Kind_DensityModel == INC_DENSITYMODEL::BOUSSINESQ) { Energy_Equation = true; @@ -4889,7 +4975,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } if (Kind_DensityModel == INC_DENSITYMODEL::VARIABLE) { - if (Kind_FluidModel != INC_IDEAL_GAS && Kind_FluidModel != INC_IDEAL_GAS_POLY && Kind_FluidModel != FLUID_MIXTURE) { + if (Kind_FluidModel != INC_IDEAL_GAS && Kind_FluidModel != INC_IDEAL_GAS_POLY && Kind_FluidModel != FLUID_MIXTURE && Kind_FluidModel != FLUID_FLAMELET) { SU2_MPI::Error("Variable density incompressible solver limited to ideal gases.\n Check the fluid model options (use INC_IDEAL_GAS, INC_IDEAL_GAS_POLY).", CURRENT_FUNCTION); } } @@ -5327,7 +5413,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i CURRENT_FUNCTION); /*--- Checks for additional species transport. ---*/ - if (Kind_Species_Model == SPECIES_MODEL::SPECIES_TRANSPORT) { + if ((Kind_Species_Model == SPECIES_MODEL::SPECIES_TRANSPORT) || (Kind_Species_Model == SPECIES_MODEL::FLAMELET)) { if (Kind_Solver != MAIN_SOLVER::INC_NAVIER_STOKES && Kind_Solver != MAIN_SOLVER::INC_RANS && Kind_Solver != MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES && @@ -5397,24 +5483,33 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i nSpecies = nSpecies_Init; /*--- Check whether some variables (or their sums) are in physical bounds. [0,1] for species related quantities. ---*/ - su2double Species_Init_Sum = 0.0; - for (unsigned short iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - checkScalarBounds(Species_Init[iSpecies], "SPECIES_INIT individual", 0.0, 1.0); - Species_Init_Sum += Species_Init[iSpecies]; - } - checkScalarBounds(Species_Init_Sum, "SPECIES_INIT sum", 0.0, 1.0); - - for (unsigned short iMarker = 0; iMarker < nMarker_Inlet_Species; iMarker++) { - su2double Inlet_SpeciesVal_Sum = 0.0; + /*--- Note, only for species transport, not for flamelet model ---*/ + if (Kind_Species_Model == SPECIES_MODEL::SPECIES_TRANSPORT) { + su2double Species_Init_Sum = 0.0; for (unsigned short iSpecies = 0; iSpecies < nSpecies; iSpecies++) { - checkScalarBounds(Inlet_SpeciesVal[iMarker][iSpecies], "MARKER_INLET_SPECIES individual", 0.0, 1.0); - Inlet_SpeciesVal_Sum += Inlet_SpeciesVal[iMarker][iSpecies]; + checkScalarBounds(Species_Init[iSpecies], "SPECIES_INIT individual", 0.0, 1.0); + Species_Init_Sum += Species_Init[iSpecies]; + } + checkScalarBounds(Species_Init_Sum, "SPECIES_INIT sum", 0.0, 1.0); + + for (unsigned short iMarker = 0; iMarker < nMarker_Inlet_Species; iMarker++) { + su2double Inlet_SpeciesVal_Sum = 0.0; + for (unsigned short iSpecies = 0; iSpecies < nSpecies; iSpecies++) { + checkScalarBounds(Inlet_SpeciesVal[iMarker][iSpecies], "MARKER_INLET_SPECIES individual", 0.0, 1.0); + Inlet_SpeciesVal_Sum += Inlet_SpeciesVal[iMarker][iSpecies]; + } + checkScalarBounds(Inlet_SpeciesVal_Sum, "MARKER_INLET_SPECIES sum", 0.0, 1.0); } - checkScalarBounds(Inlet_SpeciesVal_Sum, "MARKER_INLET_SPECIES sum", 0.0, 1.0); } } // species transport checks + /*--- Define some variables for flamelet model. ---*/ + if (Kind_Species_Model == SPECIES_MODEL::FLAMELET) { + n_control_vars = 2; + n_scalars = n_control_vars + n_user_scalars; + } + if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) { SU2_MPI::Error("BOUNDED_SCALAR discretization can only be used for incompressible problems.", CURRENT_FUNCTION); } @@ -5531,6 +5626,16 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Surface_PressureDrop = new su2double[nMarker_Analyze] (); Surface_Species_0 = new su2double[nMarker_Analyze] (); Surface_Species_Variance = new su2double[nMarker_Analyze] (); + Surface_Scalar_00 = new su2double[nMarker_Analyze] (); + Surface_Scalar_01 = new su2double[nMarker_Analyze] (); + Surface_Scalar_02 = new su2double[nMarker_Analyze] (); + Surface_Scalar_03 = new su2double[nMarker_Analyze] (); + Surface_Scalar_04 = new su2double[nMarker_Analyze] (); + Surface_Scalar_05 = new su2double[nMarker_Analyze] (); + Surface_Scalar_06 = new su2double[nMarker_Analyze] (); + Surface_Scalar_07 = new su2double[nMarker_Analyze] (); + Surface_Scalar_08 = new su2double[nMarker_Analyze] (); + Surface_Scalar_09 = new su2double[nMarker_Analyze] (); Surface_DC60 = new su2double[nMarker_Analyze] (); Surface_IDC = new su2double[nMarker_Analyze] (); Surface_IDC_Mach = new su2double[nMarker_Analyze] (); @@ -6288,6 +6393,11 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << endl; } + if (initial_PyCustom != 0) { + cout << "Using customizable initial condition" << endl; + cout << endl; + } + if (nMarker_DV != 0) { cout << "Surface(s) affected by the design variables: "; for (iMarker_DV = 0; iMarker_DV < nMarker_DV; iMarker_DV++) { @@ -8042,6 +8152,16 @@ CConfig::~CConfig() { delete[] Surface_PressureDrop; delete[] Surface_Species_0; delete[] Surface_Species_Variance; + delete[] Surface_Scalar_00; + delete[] Surface_Scalar_01; + delete[] Surface_Scalar_02; + delete[] Surface_Scalar_03; + delete[] Surface_Scalar_04; + delete[] Surface_Scalar_05; + delete[] Surface_Scalar_06; + delete[] Surface_Scalar_07; + delete[] Surface_Scalar_08; + delete[] Surface_Scalar_09; delete[] Surface_DC60; delete[] Surface_IDC; delete[] Surface_IDC_Mach; @@ -8253,6 +8373,16 @@ string CConfig::GetObjFunc_Extension(string val_filename) const { case SURFACE_PRESSURE_DROP: AdjExt = "_dp"; break; case SURFACE_SPECIES_0: AdjExt = "_avgspec0"; break; case SURFACE_SPECIES_VARIANCE: AdjExt = "_specvar"; break; + case SURFACE_SCALAR_00: AdjExt = "_avgsclr00";break; + case SURFACE_SCALAR_01: AdjExt = "_avgsclr01";break; + case SURFACE_SCALAR_02: AdjExt = "_avgsclr02";break; + case SURFACE_SCALAR_03: AdjExt = "_avgsclr03";break; + case SURFACE_SCALAR_04: AdjExt = "_avgsclr04";break; + case SURFACE_SCALAR_05: AdjExt = "_avgsclr05";break; + case SURFACE_SCALAR_06: AdjExt = "_avgsclr06";break; + case SURFACE_SCALAR_07: AdjExt = "_avgsclr07";break; + case SURFACE_SCALAR_08: AdjExt = "_avgsclr08";break; + case SURFACE_SCALAR_09: AdjExt = "_avgsclr09";break; case SURFACE_MACH: AdjExt = "_mach"; break; case CUSTOM_OBJFUNC: AdjExt = "_custom"; break; case FLOW_ANGLE_OUT: AdjExt = "_fao"; break; diff --git a/Common/src/containers/CFileReaderLUT.cpp b/Common/src/containers/CFileReaderLUT.cpp index ba3f19841eee..17f2e1ba9b61 100644 --- a/Common/src/containers/CFileReaderLUT.cpp +++ b/Common/src/containers/CFileReaderLUT.cpp @@ -367,4 +367,4 @@ bool CFileReaderLUT::GetStrippedLine(ifstream& file_stream, string& line) const /*--- return true if line is not empty, else return false ---*/ return !line.empty(); -} \ No newline at end of file +} diff --git a/Common/src/containers/CTrapezoidalMap.cpp b/Common/src/containers/CTrapezoidalMap.cpp index ef4159284379..09147d2c8b23 100644 --- a/Common/src/containers/CTrapezoidalMap.cpp +++ b/Common/src/containers/CTrapezoidalMap.cpp @@ -25,10 +25,12 @@ * License along with SU2. If not, see . */ +#include "../../Common/include/containers/CTrapezoidalMap.hpp" + #include +#include #include "../../Common/include/option_structure.hpp" -#include "../../Common/include/containers/CTrapezoidalMap.hpp" using namespace std; @@ -134,7 +136,7 @@ unsigned long CTrapezoidalMap::GetTriangle(su2double val_x, su2double val_y) { /* within that band, find edges which enclose the (val_x, val_y) point */ pair edges = GetEdges(band, val_x, val_y); - /* identify the triangle using the two edges */ + /* identify the adjacent triangles using the two edges */ std::array triangles_edge_low; for (int i = 0; i < 2; i++) triangles_edge_low[i] = edge_to_triangle[edges.first][i]; @@ -145,8 +147,7 @@ unsigned long CTrapezoidalMap::GetTriangle(su2double val_x, su2double val_y) { sort(triangles_edge_low.begin(), triangles_edge_low.end()); sort(triangles_edge_up.begin(), triangles_edge_up.end()); - // The intersection of the faces to which upper or lower belongs is - // the face that both belong to. + /* The intersection of the faces to which upper or lower belongs is the face that both belong to. */ vector triangle(1); set_intersection(triangles_edge_up.begin(), triangles_edge_up.end(), triangles_edge_low.begin(), triangles_edge_low.end(), triangle.begin()); diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 5481a8e9a977..9234019e3d80 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -441,6 +441,33 @@ class CDriver { */ unsigned long GetNumberVertices(unsigned short iMarker) const; + + /*! + * \brief Get the coordinates of the point on the mesh + * \param[in] iPoint - Point index. + * \param[in] iMesh - Mesh identifier. + * \return Number of points. + */ + vector GetCoords(unsigned long iPoint, unsigned short iMesh) const; + + /*! + * \brief Get the number of conservative state variables. + * \return Number of conservative state variables. + */ + unsigned long GetNumberStateVariables() const; + + /*! + * \brief Get the number of conservative state variables. + * \return Number of conservative state variables. + */ + unsigned long GetNumberStateVariables(const int SOLVER) const; + + /*! + * \brief Get the number of vertices in the mesh. + * \return Number of vertices. + */ + unsigned long GetNumberVertices() const; + /*! * \brief Get the number of halo vertices from a specified marker. * \param[in] iMarker - Marker identifier. @@ -462,6 +489,13 @@ class CDriver { */ unsigned long GetnTimeIter() const; + /*! + * \brief Get the number of points on the mesh + * \param[in] iMesh - Mesh identifier. + * \return Number of points. + */ + unsigned long GetnPoints(unsigned short iMesh) const; + /*! * \brief Get the current external iteration. * \return Current external iteration. @@ -488,6 +522,13 @@ class CDriver { */ unsigned long GetVertexGlobalIndex(unsigned short iMarker, unsigned long iVertex) const; + /*! + * \brief Get the global index of the solver + * \param[in] solverName - string of the solver name + * \return index of the solver + */ + unsigned long GetSolverIndex(string solverName) const; + /*! * \brief Get undeformed coordinates from the mesh solver. * \param[in] iMarker - Marker identifier. @@ -717,6 +758,34 @@ class CDriver { */ void SetHeatSource_Position(passivedouble alpha, passivedouble pos_x, passivedouble pos_y, passivedouble pos_z); + /*! + * \brief Set the conservative states at the mesh vertices. + * \param[in] values - Flow states (nPoint, nVar). + */ + void SetStates(vector> values); + + /*! + * \brief Set the conservative states at a mesh vertex. + * \param[in] iPoint - Mesh vertex index. + * \param[in] values - Flow states (nVar). + */ + void SetStates(unsigned long iPoint, vector values); + + /*! + * \brief Set the conservative states at the mesh vertices. + * \param[in] values - Flow states (nPoint, nVar). + * \param[in] SOLVER - solver type, e.g. FLOW_SOL, TURB_SOL + */ + void SetStates(vector> values, const int SOLVER); + + /*! + * \brief Set the conservative states at a mesh vertex. + * \param[in] iPoint - Mesh vertex index. + * \param[in] values - Flow states (nVar). + * \param[in] SOLVER - solver type, e.g. FLOW_SOL, TURB_SOL + */ + void SetStates(unsigned long iPoint, vector values, const int SOLVER); + /*! * \brief Set the direction of the inlet. * \param[in] iMarker - Marker index. diff --git a/SU2_CFD/include/fluid/CFluidFlamelet.hpp b/SU2_CFD/include/fluid/CFluidFlamelet.hpp new file mode 100644 index 000000000000..cf561d66baf5 --- /dev/null +++ b/SU2_CFD/include/fluid/CFluidFlamelet.hpp @@ -0,0 +1,162 @@ +/*! + * \file CFluidFlamelet.hpp + * \brief Defines the flamelet fluid model + * \author D. Mayer, T. Economon, N. Beishuizen + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "../../Common/include/containers/CLookUpTable.hpp" +#include "CFluidModel.hpp" + +class CFluidFlamelet final : public CFluidModel { + protected: + int rank; + + unsigned short n_scalars; + unsigned short n_lookups; + unsigned short n_table_sources; + unsigned short n_user_scalars; /*!< \brief number of passive reactant species. */ + unsigned short n_control_vars; /*!< \brief number of controlling variables. */ + + vector table_scalar_names; /*!< \brief vector to store names of scalar variables. */ + vector table_source_names; /*!< \brief vector to store names of scalar source variables. */ + vector table_lookup_names; /*!< \brief vector to store names of look up variables. */ + + vector table_sources; + + su2double mass_diffusivity; /*!< \brief local mass diffusivity of the mixture */ + su2double molar_weight; /*!< \brief local molar weight of the mixture */ + + vector source_scalar; + vector lookup_scalar; + + CLookUpTable* look_up_table; + + vector varnames_TD; /*!< \brief Lookup names for thermodynamic state variables. */ + vector val_vars_TD; /*!< \brief References to thermodynamic state variables. */ + + vector varnames_Sources; /*!< \brief Lookup names for scalar source terms. */ + vector val_vars_Sources; /*!< \brief References to scalar sources. */ + + vector varnames_LookUp; /*!< \brief Lookup names for passive lookup variables. */ + vector val_vars_LookUp; /*!< \brief References to lookup variables. */ + + public: + CFluidFlamelet(CConfig* config, su2double value_pressure_operating); + + ~CFluidFlamelet(); + + /*! + * \brief Set the thermodynamioc state + * \param[in] val_temperature - temperature + * \param[in] val_scalars - pointer to species mass fractions + */ + void SetTDState_T(su2double val_temperature, const su2double* val_scalars = nullptr) override; + + /*! + * \brief Set the reaction source terms for the transported species equations + * \param[in] val_scalars - pointer to species mass fractions + * \param[out] exit_code = error code + */ + unsigned long SetScalarSources(su2double* val_scalars); + + /*! + * \brief Retrieve and set the lookup values for the species + * \param[in] val_scalars - pointer to species mass fractions + */ + unsigned long SetScalarLookups(su2double* val_scalars); + + /*! + * \brief Get the total enthalpy from the tabulated temperature and species (inverse lookup) + * \param[in/out] enthalpy - total enthalpy + * \param[in] val_prog - progress variable + * \param[in] val_temp - temperature + * \param[out] exit_code = error code + */ + unsigned long GetEnthFromTemp(su2double* enthalpy, su2double val_prog, su2double val_temp, + su2double initial_value = 0); + + /*! + * \brief return a pointer to the lookup table + * \param[out] look_up_table - pointer to lookup table + */ + inline CLookUpTable* GetLookUpTable() { return look_up_table; } + + /*! + * \brief Get the mass diffusivity of the species + * \param[in] iVar - index to the species + * \param[out] mass_diffusivity - value of the mass diffusivity + */ + inline su2double GetMassDiffusivity(int iVar) final { return mass_diffusivity; } + + /*! + * \brief Get the thermal conductivity of the species + * \param[in] iVar - index to the species + * \param[out] Kt - value of the thermal conductivity + */ + inline su2double GetThermalConductivity() { return Kt; } + + /*! + * \brief Get the laminar viscosity of the species + * \param[in] iVar - index to the species + * \param[out] Mu - value of the laminar viscosity + */ + inline su2double GetLaminarViscosity() { return Mu; } + + /*! + * \brief Get the enthalpy range in the lookup table + */ + inline pair GetTableLimitsEnth() { + pair limits_Y = look_up_table->GetTableLimitsY(); + return make_pair(*limits_Y.first, *limits_Y.second); + } + + /*! + * \brief Get the progress variable range in the lookup table + */ + inline pair GetTableLimitsProg() { + pair limits_X = look_up_table->GetTableLimitsX(); + return make_pair(*limits_X.first, *limits_X.second); + } + + /*! + * \brief Get the reaction source term of a species equation + * \param[in] iVar - index to the species + */ + inline su2double GetScalarSources(int iVar) { return source_scalar[iVar]; } + + /*! + * \brief Get the reaction source term of all species equations + */ + inline su2double* GetScalarSources() { return &source_scalar[0]; } + + /*! + * \brief Get the value of the looked up variable + * \param[in] i_scalar - index to the value that we need to retrieve from the lookup table + */ + inline su2double GetScalarLookups(int i_scalar) { return lookup_scalar[i_scalar]; } + + void PreprocessLookUp(); +}; diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index ba52958f3a91..11bad7c8e297 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/SU2_CFD/include/fluid/CFluidModel.hpp @@ -32,6 +32,7 @@ #include "../../../Common/include/CConfig.hpp" #include "../../../Common/include/basic_types/datatype_structure.hpp" +#include "../../../Common/include/containers/CLookUpTable.hpp" #include "CConductivityModel.hpp" #include "CViscosityModel.hpp" #include "CDiffusivityModel.hpp" @@ -68,6 +69,7 @@ class CFluidModel { 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. */ + CLookUpTable* look_up_table; /*!< \brief the lookup table for the flamelet combustion model*/ su2double mass_diffusivity{0.0}; /*!< \brief Mass Diffusivity */ unique_ptr LaminarViscosity; /*!< \brief Laminar Viscosity Model */ @@ -137,6 +139,43 @@ class CFluidModel { */ su2double GetCv() const { return Cv; } + /*! + * \brief flamelet LUT - Get the source term of the transported scalar + */ + virtual inline su2double* GetScalarSources(){ return nullptr; } + + /*! + * \brief flamelet LUT - Get the source term of the transported scalar + * \param[in] val_ix - Index of the scalar. + */ + virtual inline su2double GetScalarSources(int val_ix){ return 0; } + + /*! + * \brief flamelet LUT - Get the number of transported scalars + */ + virtual inline unsigned short GetNScalars() {return 0; } + + /*! + * \brief flamelet LUT - Get the looked up scalar field for combustion + */ + virtual inline su2double GetScalarLookups(int){ return 0; } + + /*! + * \brief flamelet LUT - Get the actual lookup table + */ + virtual CLookUpTable* GetLookUpTable() { return look_up_table; } + + /*! + * \brief flamelet LUT - Get the total enthalpy from the temperature (reverse lookup) + */ + virtual inline unsigned long GetEnthFromTemp(su2double *enthalpy, + su2double val_prog, + su2double val_temp, + su2double initial_value=0) { return 0; } + + virtual inline pair GetTableLimitsEnth() { return make_pair(0,0); } + virtual inline pair GetTableLimitsProg() { return make_pair(0,0); } + /*! * \brief Get fluid dynamic viscosity. */ @@ -319,7 +358,14 @@ 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 Virtual member. nijso: todo: is this really necessary? + */ + virtual unsigned long SetScalarSources(su2double* val_scalars) { return 0; } + + virtual unsigned long SetScalarLookups(su2double* val_scalars) { return 0; } /*! * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity. diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index e8682cc8da58..872cec9858fb 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -368,6 +368,12 @@ class CNumerics { TransVar_j = val_transvar_j; } + /*! + * \brief Set the values of the volumetric scalar sources for the flamelet LUT approach. + * \param[in] val_scalar_sources - Values of the scalar sources. + */ + virtual void SetScalarSources(su2double *val_scalar_sources) {/* empty */} + /*! * \brief Set the gradient of the scalar variables. * \param[in] val_scalarvar_grad_i - Gradient of the scalar variable at point i. diff --git a/SU2_CFD/include/numerics/species/species_diffusion.hpp b/SU2_CFD/include/numerics/species/species_diffusion.hpp index e100496b135f..2d5279817e9d 100644 --- a/SU2_CFD/include/numerics/species/species_diffusion.hpp +++ b/SU2_CFD/include/numerics/species/species_diffusion.hpp @@ -55,6 +55,7 @@ class CAvgGrad_Species final : public CAvgGrad_Scalar { using Base::Jacobian_j; const bool turbulence; + //const bool flamelet; /*! * \brief Adds any extra variables to AD @@ -69,10 +70,12 @@ class CAvgGrad_Species final : public CAvgGrad_Scalar { * \param[in] config - Definition of the particular problem. */ void FinishResidualCalc(const CConfig* config) override { + su2double Diffusivity_Lam = 0.0; + for (auto iVar = 0u; iVar < nVar; iVar++) { + + Diffusivity_Lam = 0.5 * (Density_i * Diffusion_Coeff_i[iVar] + Density_j * Diffusion_Coeff_j[iVar]); - /* --- in case of species transport, Diffusion_Coeff is the binary diffusion coefficient --- */ - const su2double Diffusivity_Lam = 0.5 * (Density_i * Diffusion_Coeff_i[iVar] + Density_j * Diffusion_Coeff_j[iVar]); su2double Diffusivity_Turb = 0.0; if (turbulence) { diff --git a/SU2_CFD/include/numerics/species/species_sources.hpp b/SU2_CFD/include/numerics/species/species_sources.hpp index 7dd564720efb..ece7088f40fc 100644 --- a/SU2_CFD/include/numerics/species/species_sources.hpp +++ b/SU2_CFD/include/numerics/species/species_sources.hpp @@ -89,4 +89,126 @@ class CSourceAxisymmetric_Species : public CSourceBase_Species { * \return Lightweight const-view of residual and Jacobian. */ ResidualType<> ComputeResidual(const CConfig* config) override; + +}; + + +/*! + * \class CSourcePieceWise_transportedScalar_general + * \brief Class for integrating the source terms of the transported scalar turbulence model equations. + * \ingroup SourceDiscr + * \author N. Beishuizen + */ +class CSourcePieceWise_transportedScalar_general final : public CNumerics { + +private: + su2double *Residual = nullptr; + su2double **Jacobian_i = nullptr; + bool flamelet; + su2double *scalar_sources = nullptr; + + su2double Sc_t; + + su2double source_pv; + + bool incompressible; + bool viscous; + bool axisymmetric; + bool implicit; + bool inc_rans; + + /*! + * \brief Add contribution due to axisymmetric formulation to 2D residual + */ + inline void ResidualAxisymmetric(void) { + su2double yinv,Density_i,Velocity_i[3]; + if (Coord_i[1] > EPS) { + + AD::SetPreaccIn(Coord_i[1]); + + yinv = 1.0/Coord_i[1]; + + /*--- the incompressible density. Note that this is different for compressible flows ---*/ + + Density_i = V_i[nDim+2]; + + /*--- Set primitive variables at points iPoint. ---*/ + + for (auto iDim = 0u; iDim < nDim; iDim++) + Velocity_i[iDim] = V_i[iDim+1]; + + /*--- Inviscid component of the source term. ---*/ + + for (auto iVar=0u; iVar < nVar; iVar++) + Residual[iVar] -= yinv*Volume*Density_i*ScalarVar_i[iVar]*Velocity_i[1]; + + if (implicit) { + for (auto iVar=0u; iVar < nVar; iVar++) { + Jacobian_i[iVar][iVar] -= yinv*Volume*Velocity_i[1]; + } + } + + /*--- Add the viscous terms if necessary. ---*/ + + if (viscous) { + Laminar_Viscosity_i = V_i[nDim+4]; + Eddy_Viscosity_i = V_i[nDim+5]; + Thermal_Conductivity_i = V_i[nDim+6]; + + su2double Mass_Diffusivity_Tur = 0.0; + if (inc_rans) + Mass_Diffusivity_Tur = Eddy_Viscosity_i/Sc_t; + + for (auto iVar=0u; iVar < nVar; iVar++) + Residual[iVar] += yinv*Volume*(Density_i*Diffusion_Coeff_i[iVar] + Mass_Diffusivity_Tur)*ScalarVar_Grad_i[iVar][1]; + } + + } else { + + for (auto iVar=0u; iVar < nVar; iVar++) + Residual[iVar] = 0.0; + + if (implicit) { + for (auto iVar=0u; iVar < nVar; iVar++) { + for (auto jVar=0u; jVar < nVar; jVar++) + Jacobian_i[iVar][jVar] = 0.0; + } + } + + } + + } + +public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CSourcePieceWise_transportedScalar_general(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); + + /*! + * \brief Residual for source term integration. + * \param[in] config - Definition of the particular problem. + * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. + */ + ResidualType<> ComputeResidual(const CConfig* config) override; + + /*! + * \brief Set the value of the progress variable source term for the flamelet model. + * \param[in] val_sourcepv_i - Value of the source term at point i. + * \param[in] val_sourcepv_j - Value of the source term at point j. + */ + inline void SetScalarSources(su2double *val_scalar_sources) override { + for (auto i_var=0u; i_var < nVar; i_var++) + scalar_sources[i_var] = val_scalar_sources[i_var]; + } + + /*! + * \brief Destructor of the class. + */ + ~CSourcePieceWise_transportedScalar_general(void) override; + + }; diff --git a/SU2_CFD/include/output/CFlowIncOutput.hpp b/SU2_CFD/include/output/CFlowIncOutput.hpp index 7578f3a7cd27..c69e1f5650c3 100644 --- a/SU2_CFD/include/output/CFlowIncOutput.hpp +++ b/SU2_CFD/include/output/CFlowIncOutput.hpp @@ -41,6 +41,7 @@ class CFlowIncOutput final: public CFlowOutput { TURB_MODEL turb_model; /*!< \brief The kind of turbulence model*/ bool heat; /*!< \brief Boolean indicating whether have a heat problem*/ bool weakly_coupled_heat; /*!< \brief Boolean indicating whether have a weakly coupled heat equation*/ + bool flamelet; /*!< \brief Boolean indicating whether we solve the flamelet equations */ unsigned short streamwisePeriodic; /*!< \brief Boolean indicating whether it is a streamwise periodic simulation. */ bool streamwisePeriodic_temperature; /*!< \brief Boolean indicating streamwise periodic temperature is used. */ diff --git a/SU2_CFD/include/output/CFlowOutput.hpp b/SU2_CFD/include/output/CFlowOutput.hpp index 9f9b3e029b84..a7a612521ce1 100644 --- a/SU2_CFD/include/output/CFlowOutput.hpp +++ b/SU2_CFD/include/output/CFlowOutput.hpp @@ -242,6 +242,7 @@ class CFlowOutput : public CFVMOutput{ return std::stoi(std::string(s.begin() + nameLen + 1, s.end() - 1)); }; if (var.rfind("SPECIES", 0) == 0) return SPECIES_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 7); + if (var.rfind("FLAMELET", 0) == 0) return SPECIES_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 7); if (var.rfind("TURB", 0) == 0) return TURB_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 4); if (var.rfind("RAD", 0) == 0) return RAD_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 3); diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index bab8f187cae5..8295b1e837ef 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2981,6 +2981,37 @@ su2double CFVMFlowSolverBase::EvaluateCommonObjFunc(const CConfig& config) case SURFACE_SPECIES_VARIANCE: objFun += weight * config.GetSurface_Species_Variance(0); break; + case SURFACE_SCALAR_00: + objFun += weight * config.GetSurface_Scalar_00(0); + break; + case SURFACE_SCALAR_01: + objFun += weight * config.GetSurface_Scalar_01(0); + break; + case SURFACE_SCALAR_02: + objFun += weight * config.GetSurface_Scalar_02(0); + break; + case SURFACE_SCALAR_03: + objFun += weight * config.GetSurface_Scalar_03(0); + break; + case SURFACE_SCALAR_04: + objFun += weight * config.GetSurface_Scalar_04(0); + break; + case SURFACE_SCALAR_05: + objFun += weight * config.GetSurface_Scalar_05(0); + break; + case SURFACE_SCALAR_06: + objFun += weight * config.GetSurface_Scalar_06(0); + break; + case SURFACE_SCALAR_07: + objFun += weight * config.GetSurface_Scalar_07(0); + break; + case SURFACE_SCALAR_08: + objFun += weight * config.GetSurface_Scalar_08(0); + break; + case SURFACE_SCALAR_09: + objFun += weight * config.GetSurface_Scalar_09(0); + break; + case CUSTOM_OBJFUNC: objFun += weight * Total_Custom_ObjFunc; break; diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index f87fcf58c44a..e86ea9bfdfe7 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -4461,4 +4461,16 @@ class CSolver { } } + /*! + * \brief LUT premixed flamelet: Set the number of table look up misses. + * \param[in] val_n_table_miss - Number of table look up misses. + */ + inline virtual void SetNTableMisses(unsigned short val_n_table_miss) { } + + /*! + * \brief LUT premixed flamelet: Get the number of table look up misses. + * \return Number of table look up misses. + */ + inline virtual unsigned long GetNTableMisses() { return 0; } + }; diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp new file mode 100644 index 000000000000..7355030bccd2 --- /dev/null +++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp @@ -0,0 +1,204 @@ +/*! + * \file CSpeciesFlameletSolver.hpp + * \brief Headers of the CSpeciesFlameletSolver class + * \author D. Mayer, N. Beishuizen, T. Economon + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "CSpeciesSolver.hpp" + +/*! + * \class CSpeciesFlameletSolver + * \brief Main class for defining the flamelet model solver. + * \author N. Beishuizen + */ +class CSpeciesFlameletSolver final : public CSpeciesSolver { + private: + unsigned long n_table_misses; /*!< \brief number of times we failed to do a lookup from the table */ + vector conjugate_var; /*!< \brief CHT variables for each boundary and vertex. */ + + public: + /*! + * \brief Constructor. + * \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); + + /*! + * \brief Destructor of the class. + */ + ~CSpeciesFlameletSolver() = default; + + /*! + * \brief Restart residual and compute gradients. + * \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. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] Output - boolean to determine whether to print output. + */ + void Preprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output); + + /*! + * \brief Post-processing routine for the passive scalar model. + * \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); + /*! + * \brief Compute the primitive variables (diffusivities) + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] Output - boolean to determine whether to print output. + * \return - The number of non-physical points. + */ + unsigned long SetPrimitive_Variables(CSolver** solver_container, CConfig* config, bool Output); + + /*! + * \brief Set the initial condition for the scalar transport problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] ExtIter - External iteration. + */ + void SetInitialCondition(CGeometry** geometry, CSolver*** solver_container, CConfig* config, + unsigned long ExtIter) override; + + /*! + * \brief Compute the preconditioner for low-Mach flows. + * \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. + */ + void SetPreconditioner(CGeometry* geometry, CSolver** solver_container, CConfig* config); + + /*! + * \brief Source term computation. + * \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. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Source_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config, + unsigned short iMesh) override; + + /*! + * \brief Impose the Navier-Stokes wall boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_HeatFlux_Wall(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, + CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) override; + + /*! + * \brief Impose the Navier-Stokes wall boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Isothermal_Wall(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, + CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) override; + + /*! + * \brief Impose the inlet boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Inlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, + CConfig* config, unsigned short val_marker) override; + + /*! + * \brief Impose the outlet boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, + CConfig* config, unsigned short val_marker) override; + + inline void SetNTableMisses(unsigned short val_n_table_misses) override { n_table_misses = val_n_table_misses; } + + inline unsigned long GetNTableMisses() override { return n_table_misses; } + + /*! + * \brief Impose the (received) conjugate heat variables. + * \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. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_ConjugateHeat_Interface(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, CConfig* config, + unsigned short val_marker) override; + + /*! + * \brief Set the conjugate heat variables. + * \param[in] val_marker - marker index + * \param[in] val_vertex - vertex index + * \param[in] pos_var - variable position (in vector of all conjugate heat variables) + */ + inline su2double GetConjugateHeatVariable(unsigned short val_marker, unsigned long val_vertex, + unsigned short pos_var) const override { + return conjugate_var[val_marker][val_vertex][pos_var]; + } + + /*! + * \brief Set the conjugate heat variables. + * \param[in] val_marker - marker index + * \param[in] val_vertex - vertex index + * \param[in] pos_var - variable position (in vector of all conjugate heat variables) + * \param[in] relaxation factor - relaxation factor for the change of the variables + * \param[in] val_var - value of the variable + */ + inline void SetConjugateHeatVariable(unsigned short val_marker, unsigned long val_vertex, unsigned short pos_var, + su2double relaxation_factor, su2double val_var) override { + conjugate_var[val_marker][val_vertex][pos_var] = + relaxation_factor * val_var + (1.0 - relaxation_factor) * conjugate_var[val_marker][val_vertex][pos_var]; + } +}; + diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index d4d751e3ff5f..958b2918ab76 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -75,8 +75,8 @@ class CSpeciesSolver : public CScalarSolver { * \param[in] RunTime_EqSystem - System of equations which is going to be solved. * \param[in] Output - boolean to determine whether to print output. */ - void Preprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, unsigned short iMesh, - unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) final; + virtual void Preprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output); /*! * \brief Compute the viscous flux for the turbulent equation at a particular edge. @@ -100,7 +100,7 @@ class CSpeciesSolver : public CScalarSolver { * \param[in] val_marker - Surface marker where the boundary condition is applied. */ void BC_Inlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, - CConfig* config, unsigned short val_marker) final; + CConfig* config, unsigned short val_marker); /*! * \brief Store of a set of provided inlet profile values at a vertex. @@ -143,7 +143,7 @@ class CSpeciesSolver : public CScalarSolver { * \param[in] val_marker - Surface marker where the boundary condition is applied. */ void BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, - CConfig* config, unsigned short val_marker) final; + CConfig* config, unsigned short val_marker); /*--- Note that BC_Sym_Plane, BC_HeatFlux_Wall, BC_Isothermal_Wall are all treated as zero-flux BC for the * mass-factions, which can be fulfilled by no additional residual contribution on these nodes. diff --git a/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp b/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp new file mode 100644 index 000000000000..4cd73f38a69b --- /dev/null +++ b/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp @@ -0,0 +1,107 @@ +/*! + * \file CSpeciesFlameletVariable.hpp + * \brief Base class for defining the variables of the flamelet transport model. + * \author D. Mayer, T. Economon, N. Beishuizen + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "CSpeciesVariable.hpp" + +/*! + * \class CSpeciesFlameletVariable + * \brief Base class for defining the variables of the flamelet model. + */ +class CSpeciesFlameletVariable final : public CSpeciesVariable { + protected: + MatrixType source_scalar; /*!< \brief Vector of the source terms from the lookup table for each scalar equation */ + MatrixType lookup_scalar; /*!< \brief Vector of the source terms from the lookup table for each scalar equation */ + su2vector inside_table; /*!< \brief Vector of solutions inside the lookup table. */ + + public: + /*! + * \brief Constructor of the class. + * \param[in] species_inf - species variable values (initialization value). + * \param[in] npoint - Number of points/nodes/vertices in the domain. + * \param[in] ndim - Number of dimensions of the problem. + * \param[in] nvar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CSpeciesFlameletVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, + const CConfig* config); + + /*! + * \brief Set the value of the transported scalar source term + * \param[in] val_- the . + * \param[in] val_ivar - eqn. index to the . + */ + inline void SetLookupScalar(unsigned long iPoint, su2double val_lookup_scalar, unsigned short val_ivar) override { + lookup_scalar(iPoint, val_ivar) = val_lookup_scalar; + } + + /*! + + * \brief Set a source term for the specie transport equation. + * \param[in] iPoint - Node index. + * \param[in] val_ivar - Species index. + * \param[in] val_source - Source term value. + */ + inline void SetScalarSource(unsigned long iPoint, unsigned short val_ivar, su2double val_source) override { + source_scalar(iPoint, val_ivar) = val_source; + } + + /*! + * \brief Get the value of the transported scalar source term + * \param[in] val_ivar - eqn. index to the transported scalar source term + * \return Value of the progress variable source term + */ + inline su2double GetScalarSources(unsigned long iPoint, unsigned short val_ivar) const override { + return source_scalar(iPoint, val_ivar); + } + + /*! + * \brief Get the value of the looked up scalar field + * \param[in] val_ivar - eqn. index to the looked up scalar field + * \return Value of the looked up scalar field + */ + inline su2double GetScalarLookups(unsigned long iPoint, unsigned short val_ivar) const override { + return lookup_scalar(iPoint, val_ivar); + } + + /*! + * \brief Get the value of the transported scalars source term + * \return Pointer to the transported scalars source term + */ + inline su2double* GetScalarSources(unsigned long iPoint) override { return source_scalar[iPoint]; } + + /*! + * \brief Get the value of the looked up table based on the transported scalar + * \return Pointer to the transported scalars source term + */ + inline su2double* GetScalarLookups(unsigned long iPoint) override { return lookup_scalar[iPoint]; } + + inline void SetInsideTable(unsigned long iPoint, unsigned short inside) override { inside_table[iPoint] = inside; } + + inline unsigned short GetInsideTable(unsigned long iPoint) const override { return inside_table[iPoint]; } +}; diff --git a/SU2_CFD/include/variables/CSpeciesVariable.hpp b/SU2_CFD/include/variables/CSpeciesVariable.hpp index ac9b6fa44884..5932fa251cb3 100644 --- a/SU2_CFD/include/variables/CSpeciesVariable.hpp +++ b/SU2_CFD/include/variables/CSpeciesVariable.hpp @@ -38,7 +38,7 @@ class CSpeciesVariable : public CScalarVariable { MatrixType Diffusivity; /*!< \brief Matrix (nPoint,nVar) of mass diffusivities for scalar transport. */ public: - static constexpr size_t MAXNVAR = 4; /*!< \brief Max number of variables for static arrays. Increase, if necessary. */ + static constexpr size_t MAXNVAR = 20; /*!< \brief Max number of variables for static arrays. Increase, if necessary. */ /*! * \brief Constructor of the class. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index f80faebdd10f..479fa016c6bf 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2298,4 +2298,20 @@ class CVariable { virtual su2double GetSourceTerm_DispAdjoint(unsigned long iPoint, unsigned long iDim) const { return 0.0; } virtual su2double GetSourceTerm_VelAdjoint(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + /*! + * \brief LUT premixed flamelet: virtual functions for the speciesflameletvariable LUT + */ + inline virtual void SetLookupScalar(unsigned long iPoint, su2double val_lookup_scalar, unsigned short val_ivar) { } + + inline virtual void SetScalarSource(unsigned long iPoint, unsigned short val_ivar, su2double val_source) { } + + inline virtual void SetInsideTable(unsigned long iPoint, unsigned short inside) { } + + inline virtual unsigned short GetInsideTable(unsigned long iPoint) const { return 0; } + + inline virtual su2double GetScalarSources(unsigned long iPoint, unsigned short val_ivar) const { return 0.0; } + inline virtual su2double GetScalarLookups(unsigned long iPoint, unsigned short val_ivar) const { return 0.0; } + + inline virtual su2double *GetScalarSources(unsigned long iPoint) { return nullptr; } + inline virtual su2double *GetScalarLookups(unsigned long iPoint) { return nullptr; } }; diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am index 6926cadc38d7..871dc70cc834 100644 --- a/SU2_CFD/obj/Makefile.am +++ b/SU2_CFD/obj/Makefile.am @@ -52,6 +52,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/fluid/CFluidModel.cpp \ ../src/fluid/CIdealGas.cpp \ ../src/fluid/CFluidScalar.cpp \ + ../src/fluid/CFluidFlamelet.cpp \ ../src/fluid/CPengRobinson.cpp \ ../src/fluid/CVanDerWaalsGas.cpp \ ../src/fluid/CCoolProp.cpp \ @@ -180,6 +181,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/solvers/CTemplateSolver.cpp \ ../src/solvers/CTransLMSolver.cpp \ ../src/solvers/CSpeciesSolver.cpp \ + ../src/solvers/CSpeciesFlameletSolver.cpp \ ../src/solvers/CTurbSolver.cpp \ ../src/solvers/CTurbSASolver.cpp \ ../src/solvers/CTurbSSTSolver.cpp \ @@ -211,6 +213,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/variables/CIncEulerVariable.cpp \ ../src/variables/CScalarVariable.cpp \ ../src/variables/CSpeciesVariable.cpp \ + ../src/variables/CSpeciesFlameletVariable.cpp \ ../src/variables/CTurbVariable.cpp \ ../src/variables/CNSVariable.cpp \ ../src/variables/CNEMOEulerVariable.cpp \ diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 121cd94e9a55..eade62cf41c3 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1434,15 +1434,23 @@ void CDriver::InstantiateSpeciesNumerics(unsigned short nVar_Species, int offset /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + // nijso: we currently use axisymmetry with species transport using a simple inline function. if (config->GetAxisymmetric() == YES) { - numerics[iMGlevel][SPECIES_SOL][source_first_term] = new CSourceAxisymmetric_Species(nDim, nVar_Species, config); + numerics[iMGlevel][SPECIES_SOL][source_second_term] = new CSourceAxisymmetric_Species(nDim, nVar_Species, config); + } + else { + numerics[iMGlevel][SPECIES_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Species, config); + } + + if (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET){ + numerics[iMGlevel][SPECIES_SOL][source_first_term] = new CSourcePieceWise_transportedScalar_general(nDim, nVar_Species, config); } else { numerics[iMGlevel][SPECIES_SOL][source_first_term] = new CSourceNothing(nDim, nVar_Species, config); } - numerics[iMGlevel][SPECIES_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Species, config); } } + /*--- Explicit instantiation of the template above, needed because it is defined in a cpp file, instead of hpp. ---*/ template void CDriver::InstantiateSpeciesNumerics>( unsigned short, int, const CConfig*, const CSolver*, CNumerics****&) const; @@ -3031,6 +3039,8 @@ void CFluidDriver::Preprocess(unsigned long Iter) { if (config_container[iZone]->GetFluidProblem()) { for (iInst = 0; iInst < nInst[iZone]; iInst++) { solver_container[iZone][iInst][MESH_0][FLOW_SOL]->SetInitialCondition(geometry_container[iZone][INST_0], solver_container[iZone][iInst], config_container[iZone], Iter); + if (config_container[iZone]->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET) + solver_container[iZone][iInst][MESH_0][SPECIES_SOL]->SetInitialCondition(geometry_container[iZone][INST_0], solver_container[iZone][iInst], config_container[iZone], Iter); } } } diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index 99c64c233f51..697af8fd8123 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -245,6 +245,12 @@ void CMultizoneDriver::Preprocess(unsigned long TimeIter) { solver_container[iZone][INST_0], config_container[iZone], TimeIter); } + if (config_container[iZone]->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET) + solver_container[iZone][INST_0][MESH_0][SPECIES_SOL]->SetInitialCondition(geometry_container[ZONE_0][INST_0], + solver_container[ZONE_0][INST_0], + config_container[ZONE_0], TimeIter); + + } SU2_MPI::Barrier(SU2_MPI::GetComm()); diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index fe4bd854e1ab..bb0b26d47cc4 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -135,6 +135,11 @@ void CSinglezoneDriver::Preprocess(unsigned long TimeIter) { solver_container[ZONE_0][INST_0], config_container[ZONE_0], TimeIter); } + if (config_container[ZONE_0]->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET) + solver_container[ZONE_0][INST_0][MESH_0][SPECIES_SOL]->SetInitialCondition(geometry_container[ZONE_0][INST_0], + solver_container[ZONE_0][INST_0], + config_container[ZONE_0], TimeIter); + else if (config_container[ZONE_0]->GetHeatProblem()) { /*--- Set the initial condition for HEAT equation ---------------------------------------------*/ solver_container[ZONE_0][INST_0][MESH_0][HEAT_SOL]->SetInitialCondition(geometry_container[ZONE_0][INST_0], diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp new file mode 100644 index 000000000000..a3426cca8cae --- /dev/null +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -0,0 +1,231 @@ +/*! + * \file CfluidFlamelet.cpp + * \brief Main subroutines of CFluidFlamelet class + * \author D. Mayer, T. Economon, N. Beishuizen + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../include/fluid/CFluidFlamelet.hpp" + +#include "../../../Common/include/containers/CLookUpTable.hpp" + +CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operating) : CFluidModel() { +#ifdef HAVE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &rank); +#endif + + /* -- number of auxiliary species transport equations: 1=CO, 2=NOx --- */ + n_user_scalars = config->GetNUserScalars(); + n_control_vars = config->GetNControlVars(); + n_scalars = config->GetNScalars(); + + if (rank == MASTER_NODE) { + cout << "Number of scalars: " << n_scalars << endl; + cout << "Number of user scalars: " << n_user_scalars << endl; + cout << "Number of control variables: " << n_control_vars << endl; + } + + if (rank == MASTER_NODE) { + cout << "*****************************************" << endl; + cout << "*** initializing the lookup table ***" << endl; + cout << "*****************************************" << endl; + } + + table_scalar_names.resize(n_scalars); + table_scalar_names[I_ENTH] = "EnthalpyTot"; + table_scalar_names[I_PROGVAR] = "ProgressVariable"; + /*--- auxiliary species transport equations---*/ + for (size_t i_aux = 0; i_aux < n_user_scalars; i_aux++) { + table_scalar_names[n_control_vars + i_aux] = config->GetUserScalarName(i_aux); + } + + config->SetLUTScalarNames(table_scalar_names); + + /*--- we currently only need 1 source term from the LUT for the progress variable + and each auxiliary equations needs 2 source terms ---*/ + n_table_sources = 1 + 2 * n_user_scalars; + config->SetNLUTSources(n_table_sources); + + table_source_names.resize(n_table_sources); + table_sources.resize(n_table_sources); + table_source_names[I_SRC_TOT_PROGVAR] = "ProdRateTot-PV"; + /*--- No source term for enthalpy ---*/ + + /*--- For the auxiliary equations, we use a positive (production) and a negative (consumption) term: + S_tot = S_PROD + S_CONS * Y ---*/ + + for (size_t i_aux = 0; i_aux < n_user_scalars; i_aux++) { + /*--- Order of the source terms: S_prod_1, S_cons_1, S_prod_2, S_cons_2, ...---*/ + table_source_names[1 + 2 * i_aux] = config->GetUserSourceName(2 * i_aux); + table_source_names[1 + 2 * i_aux + 1] = config->GetUserSourceName(2 * i_aux + 1); + } + + config->SetLUTSourceNames(table_source_names); + + look_up_table = new CLookUpTable(config->GetFileNameLUT(), table_scalar_names[I_PROGVAR], table_scalar_names[I_ENTH]); + + n_lookups = config->GetNLookups(); + table_lookup_names.resize(n_lookups); + for (int i_lookup = 0; i_lookup < n_lookups; ++i_lookup) { + table_lookup_names[i_lookup] = config->GetLUTLookupName(i_lookup); + } + + source_scalar.resize(n_scalars); + lookup_scalar.resize(n_lookups); + + Pressure = value_pressure_operating; + + PreprocessLookUp(); +} + +CFluidFlamelet::~CFluidFlamelet() { + if (look_up_table != NULL) delete look_up_table; +} + +/*--- do a lookup for the list of variables in table_lookup_names, for visualization purposes ---*/ +unsigned long CFluidFlamelet::SetScalarLookups(su2double* val_scalars) { + su2double enth = val_scalars[I_ENTH]; + su2double prog = val_scalars[I_PROGVAR]; + + /*--- perform table look ups ---*/ + unsigned long exit_code = look_up_table->LookUp_XY(table_lookup_names, lookup_scalar, prog, enth); + + return exit_code; +} + +/*--- set the source terms for the transport equations ---*/ +unsigned long CFluidFlamelet::SetScalarSources(su2double* val_scalars) { + table_sources[0] = 0.0; + + /*--- value for the progress variable and enthalpy ---*/ + su2double enth = val_scalars[I_ENTH]; + su2double prog = val_scalars[I_PROGVAR]; + + /*--- perform table look ups ---*/ + unsigned long exit_code = look_up_table->LookUp_XY(varnames_Sources, val_vars_Sources, prog, enth); + + /*--- the source term for the progress variable is always positive by construction, but we clip from below --- */ + source_scalar[I_PROGVAR] = max(EPS, table_sources[I_SRC_TOT_PROGVAR]); + source_scalar[I_ENTH] = 0.0; + + /*--- source term for the auxiliary species transport equations---*/ + for (size_t i_aux = 0; i_aux < n_user_scalars; i_aux++) { + /*--- The source term for the auxiliary equations consists of a production term and a consumption term: + S_TOT = S_PROD + S_CONS * Y ---*/ + su2double y_aux = val_scalars[n_control_vars + i_aux]; + su2double source_prod = table_sources[1 + 2 * i_aux]; + su2double source_cons = table_sources[1 + 2 * i_aux + 1]; + source_scalar[n_control_vars + i_aux] = source_prod + source_cons * y_aux; + } + + return exit_code; +} + +void CFluidFlamelet::SetTDState_T(su2double val_temperature, const su2double* val_scalars) { + su2double val_enth = val_scalars[I_ENTH]; + su2double val_prog = val_scalars[I_PROGVAR]; + + /*--- add all quantities and their address to the look up vectors ---*/ + look_up_table->LookUp_XY(varnames_TD, val_vars_TD, val_prog, val_enth); + + /*--- compute Cv from Cp and molar weight of the mixture (ideal gas) ---*/ + Cv = Cp - UNIVERSAL_GAS_CONSTANT / molar_weight; +} + +unsigned long CFluidFlamelet::GetEnthFromTemp(su2double* val_enth, su2double val_prog, su2double val_temp, + su2double initial_value) { + su2double delta_temp_final = 0.01; /* convergence criterion for temperature in [K] */ + su2double enth_iter = initial_value; + su2double delta_enth; + su2double delta_temp_iter = 1e10; + unsigned long exit_code = 0; + vector look_up_tags; + vector look_up_data; + int counter_limit = 50; + + int counter = 0; + while ((abs(delta_temp_iter) > delta_temp_final) && (counter++ < counter_limit)) { + /* look up temperature and heat capacity */ + look_up_table->LookUp_XY(varnames_TD, val_vars_TD, val_prog, enth_iter); + + /*--- calculate delta_temperature ---*/ + delta_temp_iter = val_temp - Temperature; + + /* calculate delta_enthalpy following dh = cp * dT */ + delta_enth = Cp * delta_temp_iter; + + /*--- update enthalpy ---*/ + enth_iter += delta_enth; + } + + /*--- set enthalpy value ---*/ + *val_enth = enth_iter; + + if (counter >= counter_limit) { + exit_code = 1; + } + + return exit_code; +} + +void CFluidFlamelet::PreprocessLookUp() { + /*--- Set lookup names and variables for all relevant lookup processes in the fluid model ---*/ + + /*--- Thermodynamic state variables ---*/ + varnames_TD.resize(7); + val_vars_TD.resize(7); + + /*--- The string in varnames_TD is the actual string as it appears in the LUT file ---*/ + varnames_TD[0] = "Temperature"; + val_vars_TD[0] = &Temperature; + varnames_TD[1] = "Density"; + val_vars_TD[1] = &Density; + varnames_TD[2] = "Cp"; + val_vars_TD[2] = &Cp; + varnames_TD[3] = "ViscosityDyn"; + val_vars_TD[3] = Μ + varnames_TD[4] = "Conductivity"; + val_vars_TD[4] = &Kt; + varnames_TD[5] = "DiffusionCoefficient"; + val_vars_TD[5] = &mass_diffusivity; + varnames_TD[6] = "MolarWeightMix"; + val_vars_TD[6] = &molar_weight; + + /*--- Source term variables ---*/ + varnames_Sources.resize(n_table_sources); + val_vars_Sources.resize(n_table_sources); + + for (size_t iSource = 0; iSource < n_table_sources; iSource++) { + varnames_Sources[iSource] = table_source_names[iSource]; + val_vars_Sources[iSource] = &table_sources[iSource]; + } + + /*--- Passive lookups ---*/ + varnames_LookUp.resize(n_lookups); + val_vars_LookUp.resize(n_lookups); + + for (size_t iLookup = 0; iLookup < n_lookups; iLookup++) { + varnames_LookUp[iLookup] = table_lookup_names[iLookup]; + val_vars_LookUp[iLookup] = &lookup_scalar[iLookup]; + } +} diff --git a/SU2_CFD/src/fluid/CFluidModel.cpp b/SU2_CFD/src/fluid/CFluidModel.cpp index 15ab39b24e22..9a14ab653f83 100644 --- a/SU2_CFD/src/fluid/CFluidModel.cpp +++ b/SU2_CFD/src/fluid/CFluidModel.cpp @@ -41,6 +41,7 @@ #include "../../include/fluid/CPolynomialConductivityRANS.hpp" #include "../../include/fluid/CPolynomialViscosity.hpp" #include "../../include/fluid/CSutherland.hpp" +#include "../../include/fluid/CFluidFlamelet.hpp" #include "../../include/fluid/CCoolPropViscosity.hpp" #include "../../include/fluid/CConstantLewisDiffusivity.hpp" #include "../../include/fluid/CCoolPropConductivity.hpp" @@ -58,6 +59,9 @@ unique_ptr CFluidModel::MakeLaminarViscosityModel(const CConfig case VISCOSITYMODEL::POLYNOMIAL: return unique_ptr>( new CPolynomialViscosity(config->GetMu_PolyCoeffND())); + case VISCOSITYMODEL::FLAMELET: + /*--- Viscosity is obtained from the LUT ---*/ + return nullptr; default: SU2_MPI::Error("Viscosity model not available.", CURRENT_FUNCTION); return nullptr; @@ -108,10 +112,12 @@ unique_ptr CFluidModel::MakeThermalConductivityModel(const C new CPolynomialConductivity(config->GetKt_PolyCoeffND())); } break; + case CONDUCTIVITYMODEL::FLAMELET: + /*--- conductivity is obtained from the LUT ---*/ + return nullptr; default: SU2_MPI::Error("Conductivity model not available.", CURRENT_FUNCTION); return nullptr; - break; } } @@ -134,6 +140,10 @@ unique_ptr CFluidModel::MakeMassDiffusivityModel(const CConfi return unique_ptr( new CConstantLewisDiffusivity(config->GetConstant_Lewis_Number(iSpecies))); break; + case DIFFUSIVITYMODEL::FLAMELET: + /* do nothing. Diffusivity is obtained from the table and set in setTDState_T */ + return nullptr; + break; default: SU2_MPI::Error("Diffusivity model not available.", CURRENT_FUNCTION); return nullptr; diff --git a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp index 56fcea263e20..594863bb9a1d 100644 --- a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp +++ b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp @@ -104,6 +104,7 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet switch (donor_config->GetKind_ConductivityModel()) { case CONDUCTIVITYMODEL::CONSTANT: + case CONDUCTIVITYMODEL::FLAMELET: case CONDUCTIVITYMODEL::COOLPROP: thermal_conductivity = thermal_conductivityND*donor_config->GetThermal_Conductivity_Ref(); break; diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index fd5a5252d0e0..1f190e5a23a5 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -6,6 +6,7 @@ su2_cfd_src = files(['definition_structure.cpp', su2_cfd_src += files(['fluid/CFluidModel.cpp', 'fluid/CIdealGas.cpp', 'fluid/CFluidScalar.cpp', + 'fluid/CFluidFlamelet.cpp', 'fluid/CPengRobinson.cpp', 'fluid/CVanDerWaalsGas.cpp', 'fluid/CCoolProp.cpp', @@ -60,6 +61,7 @@ su2_cfd_src += files(['variables/CIncNSVariable.cpp', 'variables/CTurbVariable.cpp', 'variables/CScalarVariable.cpp', 'variables/CSpeciesVariable.cpp', + 'variables/CSpeciesFlameletVariable.cpp', 'variables/CAdjNSVariable.cpp', 'variables/CBaselineVariable.cpp', 'variables/CDiscAdjFEABoundVariable.cpp', @@ -109,6 +111,7 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp', 'solvers/CSolver.cpp', 'solvers/CTemplateSolver.cpp', 'solvers/CSpeciesSolver.cpp', + 'solvers/CSpeciesFlameletSolver.cpp', 'solvers/CTransLMSolver.cpp', 'solvers/CTurbSolver.cpp', 'solvers/CTurbSASolver.cpp', diff --git a/SU2_CFD/src/numerics/species/species_sources.cpp b/SU2_CFD/src/numerics/species/species_sources.cpp index 06a26902b8a7..b435fa66e15c 100644 --- a/SU2_CFD/src/numerics/species/species_sources.cpp +++ b/SU2_CFD/src/numerics/species/species_sources.cpp @@ -140,3 +140,75 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Species::ComputeResidual(const template class CSourceAxisymmetric_Species >; template class CSourceAxisymmetric_Species >; template class CSourceAxisymmetric_Species >; + + +CSourcePieceWise_transportedScalar_general::CSourcePieceWise_transportedScalar_general(unsigned short val_nDim, + unsigned short val_nVar, + const CConfig* config) : + CNumerics(val_nDim, val_nVar, config) { + + incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); + axisymmetric = config->GetAxisymmetric(); + viscous = config->GetViscous(); + implicit = (config->GetKind_TimeIntScheme_Species() == EULER_IMPLICIT); + flamelet = (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET); + inc_rans = (config->GetKind_Solver() == MAIN_SOLVER::INC_RANS) || (config->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_INC_RANS); + + Sc_t = config->GetSchmidt_Number_Turbulent(); + + Residual = new su2double [nVar]; + scalar_sources = new su2double [nVar]; + Jacobian_i = new su2double* [nVar]; + + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double [nVar] (); + } +} + +CSourcePieceWise_transportedScalar_general::~CSourcePieceWise_transportedScalar_general(void){ + delete [] Residual; + delete [] scalar_sources; + if (Jacobian_i != nullptr) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + delete [] Jacobian_i[iVar]; + } + delete [] Jacobian_i; + } +} + +CNumerics::ResidualType<> CSourcePieceWise_transportedScalar_general::ComputeResidual(const CConfig* config) { + + AD::StartPreacc(); + AD::SetPreaccIn(ScalarVar_i, nVar); + AD::SetPreaccIn(scalar_sources, nVar); + AD::SetPreaccIn(Volume); + + if (incompressible) { + AD::SetPreaccIn(V_i, nDim+6); + + } + else { + AD::SetPreaccIn(V_i, nDim+7); + + } + + /*--- Implicit part for production term (to do). ---*/ + for (auto i_var = 0; i_var < nVar; i_var++) { + Residual[i_var] = scalar_sources[i_var] * Volume; + for (auto j_var = 0; j_var < nVar; j_var++) { + Jacobian_i[i_var][j_var] = 0.0; + } + } + + /*--- Contribution due to 2D axisymmetric formulation ---*/ + + if (axisymmetric) ResidualAxisymmetric(); + + /*--- Implicit part ---*/ + + AD::SetPreaccOut(Residual, nVar); + AD::EndPreacc(); + + return ResidualType<>(Residual, Jacobian_i, nullptr); + +} diff --git a/SU2_CFD/src/output/CAdjFlowIncOutput.cpp b/SU2_CFD/src/output/CAdjFlowIncOutput.cpp index 74dfbebcfd1f..002551ed9658 100644 --- a/SU2_CFD/src/output/CAdjFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CAdjFlowIncOutput.cpp @@ -109,9 +109,9 @@ void CAdjFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: Root-mean square residual of the adjoint Velocity y-component. AddHistoryOutput("RMS_ADJ_VELOCITY-Y", "rms[A_V]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint Velocity y-component.", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Root-mean square residual of the adjoint Velocity z-component. - AddHistoryOutput("RMS_ADJ_VELOCITY-Z", "rms[A_W]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint Velocity z-component.", HistoryFieldType::RESIDUAL); + if (nDim == 3) AddHistoryOutput("RMS_ADJ_VELOCITY-Z", "rms[A_W]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint Velocity z-component.", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the temperature. - AddHistoryOutput("RMS_ADJ_TEMPERATURE", "rms[A_T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + if (heat || weakly_coupled_heat) AddHistoryOutput("RMS_ADJ_TEMPERATURE", "rms[A_T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { AddHistoryOutput("ADJOINT_SOLEXTRA", "Adjoint_SolExtra", ScreenOutputFormat::FIXED, "ADJOINT_SOLEXTRA", "Adjoint value of the first extra Solution.", HistoryFieldType::COEFFICIENT); @@ -133,9 +133,9 @@ void CAdjFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: Maximum residual of the adjoint Velocity y-component AddHistoryOutput("MAX_ADJ_VELOCITY-Y", "max[A_RhoV]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint Velocity y-component", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the adjoint Velocity z-component - AddHistoryOutput("MAX_ADJ_VELOCITY-Z", "max[A_RhoW]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint Velocity z-component", HistoryFieldType::RESIDUAL); + if (nDim == 3) AddHistoryOutput("MAX_ADJ_VELOCITY-Z", "max[A_RhoW]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint Velocity z-component", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the temperature. - AddHistoryOutput("MAX_ADJ_TEMPERATURE", "max[A_T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature.", HistoryFieldType::RESIDUAL); + if (heat || weakly_coupled_heat) AddHistoryOutput("MAX_ADJ_TEMPERATURE", "max[A_T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature.", HistoryFieldType::RESIDUAL); AddHistoryOutputFields_AdjScalarMAX_RES(config); /// END_GROUP @@ -148,9 +148,9 @@ void CAdjFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: BGS residual of the adjoint Velocity y-component AddHistoryOutput("BGS_ADJ_VELOCITY-Y", "bgs[A_RhoV]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Velocity y-component", HistoryFieldType::RESIDUAL); /// DESCRIPTION: BGS residual of the adjoint Velocity z-component - AddHistoryOutput("BGS_ADJ_VELOCITY-Z", "bgs[A_RhoW]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Velocity z-component", HistoryFieldType::RESIDUAL); + if (nDim == 3) AddHistoryOutput("BGS_ADJ_VELOCITY-Z", "bgs[A_RhoW]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Velocity z-component", HistoryFieldType::RESIDUAL); /// DESCRIPTION: BGS residual of the temperature. - AddHistoryOutput("BGS_ADJ_TEMPERATURE", "bgs[A_T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + if (heat || weakly_coupled_heat) AddHistoryOutput("BGS_ADJ_TEMPERATURE", "bgs[A_T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); AddHistoryOutputFields_AdjScalarBGS_RES(config); diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index ec2b2e108d2e..b2ac146548e7 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -38,6 +38,7 @@ CFlowIncOutput::CFlowIncOutput(CConfig *config, unsigned short nDim) : CFlowOutp heat = config->GetEnergy_Equation(); weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); + flamelet = (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET); streamwisePeriodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE); streamwisePeriodic_temperature = config->GetStreamwise_Periodic_Temperature(); @@ -295,7 +296,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("VELOCITY-Y", "Velocity_y", "SOLUTION", "y-component of the velocity vector"); if (nDim == 3) AddVolumeOutput("VELOCITY-Z", "Velocity_z", "SOLUTION", "z-component of the velocity vector"); - if (heat || weakly_coupled_heat) + if (heat || weakly_coupled_heat || flamelet) AddVolumeOutput("TEMPERATURE", "Temperature","SOLUTION", "Temperature"); SetVolumeOutputFields_ScalarSolution(config); @@ -350,7 +351,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("LIMITER_VELOCITY-Y", "Limiter_Velocity_y", "LIMITER", "Limiter value of the y-velocity"); if (nDim == 3) AddVolumeOutput("LIMITER_VELOCITY-Z", "Limiter_Velocity_z", "LIMITER", "Limiter value of the z-velocity"); - if (heat || weakly_coupled_heat) + if (heat || weakly_coupled_heat || flamelet) AddVolumeOutput("LIMITER_TEMPERATURE", "Limiter_Temperature", "LIMITER", "Limiter value of the temperature"); } @@ -385,7 +386,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve if (nDim == 3) SetVolumeOutputValue("VELOCITY-Z", iPoint, Node_Flow->GetSolution(iPoint, 3)); - if (heat) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetSolution(iPoint, nDim+1)); + if (heat || flamelet) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetSolution(iPoint, nDim+1)); if (weakly_coupled_heat) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0)); // Radiation solver diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index d24bedd03bf9..44877ffb6499 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -27,6 +27,8 @@ #include #include +#include +#include #include "../../include/output/CFlowOutput.hpp" @@ -86,6 +88,16 @@ void CFlowOutput::AddAnalyzeSurfaceOutput(const CConfig *config){ } /// DESCRIPTION: Species Variance AddHistoryOutput("SURFACE_SPECIES_VARIANCE", "Species_Variance", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF", "Total species variance, measure for mixing quality. On all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); + + if (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET) { + /// DESCRIPTION: Average flamelet user scalars + for (unsigned short i_var = 0; i_var < config->GetNScalars(); i_var++) { + std::stringstream str_i_var; + str_i_var << std::setw(2) << std::setfill('0') << i_var; + AddHistoryOutput("SURFACE_SCALAR_" + str_i_var.str(), "Avg_Scalar_" + str_i_var.str(), ScreenOutputFormat::FIXED, "FLAMELET_COEFF_SURF", "Average of scalar " + std::to_string(i_var) + " on all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); + } + } + } /// END_GROUP @@ -128,6 +140,16 @@ void CFlowOutput::AddAnalyzeSurfaceOutput(const CConfig *config){ } /// DESCRIPTION: Species Variance AddHistoryOutputPerSurface("SURFACE_SPECIES_VARIANCE", "Species_Variance", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); + + if (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET) { + /// DESCRIPTION: Average flamelet user scalars + for (unsigned short i_var = 0; i_var < config->GetNScalars(); i_var++) { + std::stringstream str_i_var; + str_i_var << std::setw(2) << std::setfill('0') << i_var; + AddHistoryOutputPerSurface("SURFACE_SCALAR_" + str_i_var.str(), "Avg_Scalar_" + str_i_var.str(), ScreenOutputFormat::FIXED, "FLAMELET_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); + } + } + } /// END_GROUP } @@ -152,13 +174,16 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry const bool energy = config->GetEnergy_Equation(); const bool streamwisePeriodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE); const bool species = config->GetKind_Species_Model() != SPECIES_MODEL::NONE; + const bool flamelet = config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET; const auto nSpecies = config->GetnSpecies(); + const auto nScalars = config->GetNScalars(); const bool axisymmetric = config->GetAxisymmetric(); const unsigned short nMarker_Analyze = config->GetnMarker_Analyze(); const auto flow_nodes = solver[FLOW_SOL]->GetNodes(); const CVariable* species_nodes = species ? solver[SPECIES_SOL]->GetNodes() : nullptr; + const CVariable* scalar_nodes = flamelet ? solver[SPECIES_SOL]->GetNodes() : nullptr; vector Surface_MassFlow (nMarker,0.0); vector Surface_Mach (nMarker,0.0); @@ -176,6 +201,8 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry vector Surface_MassFlow_Abs (nMarker,0.0); su2activematrix Surface_Species(nMarker, nSpecies); Surface_Species = su2double(0.0); + su2activematrix Surface_Scalars(nMarker, nScalars); + Surface_Scalars = su2double(0.0); su2double Tot_Surface_MassFlow = 0.0; su2double Tot_Surface_Mach = 0.0; @@ -191,6 +218,7 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry su2double Tot_Momentum_Distortion = 0.0; su2double Tot_SecondOverUniformity = 0.0; vector Tot_Surface_Species(nSpecies,0.0); + vector Tot_Surface_Scalar(nScalars,0.0); /*--- Compute the numerical fan face Mach number, and the total area of the inflow ---*/ @@ -275,6 +303,10 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry for (unsigned short iVar = 0; iVar < nSpecies; iVar++) Surface_Species(iMarker, iVar) += species_nodes->GetSolution(iPoint, iVar)*Weight; + if (flamelet) + for (unsigned short iVar = 0; iVar < config->GetNScalars(); iVar++) + Surface_Scalars(iMarker, iVar) += scalar_nodes->GetSolution(iPoint, iVar)*Weight; + /*--- For now, always used the area to weight the uniformities. ---*/ Weight = abs(Area); @@ -304,6 +336,8 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry vector Surface_MassFlow_Abs_Local (nMarker_Analyze,0.0); su2activematrix Surface_Species_Local(nMarker_Analyze,nSpecies); Surface_Species_Local = su2double(0.0); + su2activematrix Surface_Scalars_Local(nMarker_Analyze,nScalars); + Surface_Scalars_Local = su2double(0.0); vector Surface_MassFlow_Total (nMarker_Analyze,0.0); vector Surface_Mach_Total (nMarker_Analyze,0.0); @@ -320,6 +354,8 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry vector Surface_MassFlow_Abs_Total (nMarker_Analyze,0.0); su2activematrix Surface_Species_Total(nMarker_Analyze,nSpecies); Surface_Species_Total = su2double(0.0); + su2activematrix Surface_Scalars_Total(nMarker_Analyze,nScalars); + Surface_Scalars_Total = su2double(0.0); vector Surface_MomentumDistortion_Total (nMarker_Analyze,0.0); @@ -349,6 +385,8 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry Surface_MassFlow_Abs_Local[iMarker_Analyze] += Surface_MassFlow_Abs[iMarker]; for (unsigned short iVar = 0; iVar < nSpecies; iVar++) Surface_Species_Local(iMarker_Analyze, iVar) += Surface_Species(iMarker, iVar); + for (unsigned short iVar = 0; iVar < nScalars; iVar++) + Surface_Scalars_Local(iMarker_Analyze, iVar) += Surface_Scalars(iMarker, iVar); } } @@ -379,7 +417,7 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry Allreduce(Surface_Area_Local, Surface_Area_Total); Allreduce(Surface_MassFlow_Abs_Local, Surface_MassFlow_Abs_Total); Allreduce_su2activematrix(Surface_Species_Local, Surface_Species_Total); - + Allreduce_su2activematrix(Surface_Scalars_Local, Surface_Scalars_Total); /*--- Compute the value of Surface_Area_Total, and Surface_Pressure_Total, and set the value in the config structure for future use ---*/ @@ -401,6 +439,8 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry Surface_TotalPressure_Total[iMarker_Analyze] /= Weight; for (unsigned short iVar = 0; iVar < nSpecies; iVar++) Surface_Species_Total(iMarker_Analyze, iVar) /= Weight; + for (unsigned short iVar = 0; iVar < nScalars; iVar++) + Surface_Scalars_Total(iMarker_Analyze, iVar) /= Weight; } else { Surface_Mach_Total[iMarker_Analyze] = 0.0; @@ -413,6 +453,8 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry Surface_TotalPressure_Total[iMarker_Analyze] = 0.0; for (unsigned short iVar = 0; iVar < nSpecies; iVar++) Surface_Species_Total(iMarker_Analyze, iVar) = 0.0; + for (unsigned short iVar = 0; iVar < nScalars; iVar++) + Surface_Scalars_Total(iMarker_Analyze, iVar) = 0.0; } /*--- Compute flow uniformity parameters separately (always area for now). ---*/ @@ -512,6 +554,36 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry config->SetSurface_Species_0(iMarker_Analyze, Species); } } + + if (flamelet) { + for (unsigned short i_var = 0; i_var < nScalars; i_var++) { + su2double scalar = Surface_Scalars_Total(iMarker_Analyze, i_var); + std::stringstream str_i_var; + str_i_var << std::setw(2) << std::setfill('0') << i_var; + SetHistoryOutputPerSurfaceValue("SURFACE_SCALAR_" + str_i_var.str(), scalar, iMarker_Analyze); + Tot_Surface_Scalar[i_var] += scalar; + if (i_var == 0) + config->SetSurface_Scalar_00(iMarker_Analyze, scalar); + if (i_var == 1) + config->SetSurface_Scalar_01(iMarker_Analyze, scalar); + if (i_var == 2) + config->SetSurface_Scalar_02(iMarker_Analyze, scalar); + if (i_var == 3) + config->SetSurface_Scalar_03(iMarker_Analyze, scalar); + if (i_var == 4) + config->SetSurface_Scalar_04(iMarker_Analyze, scalar); + if (i_var == 5) + config->SetSurface_Scalar_05(iMarker_Analyze, scalar); + if (i_var == 6) + config->SetSurface_Scalar_06(iMarker_Analyze, scalar); + if (i_var == 7) + config->SetSurface_Scalar_07(iMarker_Analyze, scalar); + if (i_var == 8) + config->SetSurface_Scalar_08(iMarker_Analyze, scalar); + if (i_var == 9) + config->SetSurface_Scalar_09(iMarker_Analyze, scalar); + } + } } /*--- Compute the average static pressure drop between two surfaces. Note @@ -549,6 +621,14 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry Surface_Area_Total); } + if (flamelet) { + for (unsigned short i_var = 0; i_var < nScalars; i_var++){ + std::stringstream str_i_var; + str_i_var << std::setw(2) << std::setfill('0') << i_var; + SetHistoryOutputValue("SURFACE_SCALAR_" + str_i_var.str(), Tot_Surface_Scalar[i_var]); + } + } + if ((rank == MASTER_NODE) && !config->GetDiscrete_Adjoint() && output) { cout.precision(6); @@ -922,10 +1002,24 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) { case TURB_TRANS_MODEL::NONE: break; } - if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { - for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - AddHistoryOutput("RMS_SPECIES_" + std::to_string(iVar), "rms[rho*Y_" + std::to_string(iVar)+"]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of transported species.", HistoryFieldType::RESIDUAL); + switch (config->GetKind_Species_Model()) { + case SPECIES_MODEL::SPECIES_TRANSPORT: { + for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { + AddHistoryOutput("RMS_SPECIES_" + std::to_string(iVar), "rms[rho*Y_" + std::to_string(iVar)+"]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of transported species.", HistoryFieldType::RESIDUAL); + } + break; + } + case SPECIES_MODEL::FLAMELET: { + AddHistoryOutput("RMS_PROGRESS_VARIABLE", "rms[PV]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the progress variable equation.", HistoryFieldType::RESIDUAL); + AddHistoryOutput("RMS_TOTAL_ENTHALPY", "rms[Enth]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the total enthalpy equation.", HistoryFieldType::RESIDUAL); + /*--- auxiliary species transport ---*/ + for(auto i_scalar=0u; i_scalar < config->GetNUserScalars(); i_scalar++){ + string scalar_name = config->GetUserScalarName(i_scalar); + AddHistoryOutput("RMS_"+scalar_name, "rms["+scalar_name+"]", ScreenOutputFormat::FIXED , "RMS_RES", "Root-mean squared residual of the "+scalar_name+" mass fraction equation." , HistoryFieldType::RESIDUAL); + } + break; } + case SPECIES_MODEL::NONE: break; } } @@ -1068,18 +1162,37 @@ void CFlowOutput::LoadHistoryData_Scalar(const CConfig* config, const CSolver* c case TURB_TRANS_MODEL::NONE: break; } - if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { - for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - SetHistoryOutputValue("RMS_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_RMS(iVar))); - SetHistoryOutputValue("MAX_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_Max(iVar))); - if (multiZone) { - SetHistoryOutputValue("BGS_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_BGS(iVar))); + switch(config->GetKind_Species_Model()) { + case SPECIES_MODEL::SPECIES_TRANSPORT: { + for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { + SetHistoryOutputValue("RMS_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_RMS(iVar))); + SetHistoryOutputValue("MAX_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_Max(iVar))); + if (multiZone) { + SetHistoryOutputValue("BGS_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_BGS(iVar))); + } + } + SetHistoryOutputValue("LINSOL_ITER_SPECIES", solver[SPECIES_SOL]->GetIterLinSolver()); + SetHistoryOutputValue("LINSOL_RESIDUAL_SPECIES", log10(solver[SPECIES_SOL]->GetResLinSolver())); + break; + } + + case SPECIES_MODEL::FLAMELET: { + SetHistoryOutputValue("RMS_PROGRESS_VARIABLE", log10(solver[SPECIES_SOL]->GetRes_RMS(I_PROGVAR))); + SetHistoryOutputValue("RMS_TOTAL_ENTHALPY", log10(solver[SPECIES_SOL]->GetRes_RMS(I_ENTH))); + /*--- auxiliary species transport ---*/ + for (unsigned short iReactant=0; iReactantGetNUserScalars(); iReactant++){ + string species_name = config->GetUserScalarName(iReactant); + SetHistoryOutputValue("RMS_" + species_name, log10(solver[SPECIES_SOL]->GetRes_RMS(config->GetNControlVars() + iReactant))); } + + SetHistoryOutputValue("LINSOL_ITER_SPECIES", solver[SPECIES_SOL]->GetIterLinSolver()); + SetHistoryOutputValue("LINSOL_RESIDUAL_SPECIES", log10(solver[SPECIES_SOL]->GetResLinSolver())); + break; } - SetHistoryOutputValue("LINSOL_ITER_SPECIES", solver[SPECIES_SOL]->GetIterLinSolver()); - SetHistoryOutputValue("LINSOL_RESIDUAL_SPECIES", log10(solver[SPECIES_SOL]->GetResLinSolver())); + case SPECIES_MODEL::NONE: break; } + } void CFlowOutput::SetVolumeOutputFields_ScalarSolution(const CConfig* config){ @@ -1110,9 +1223,39 @@ void CFlowOutput::SetVolumeOutputFields_ScalarSolution(const CConfig* config){ break; } - if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { - for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) - AddVolumeOutput("SPECIES_" + std::to_string(iVar), "Species_" + std::to_string(iVar), "SOLUTION", "Species_" + std::to_string(iVar) + " mass fraction"); + switch (config->GetKind_Species_Model()) { + case SPECIES_MODEL::SPECIES_TRANSPORT: + for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) + AddVolumeOutput("SPECIES_" + std::to_string(iVar), "Species_" + std::to_string(iVar), "SOLUTION", "Species_" + std::to_string(iVar) + " mass fraction"); + break; + case SPECIES_MODEL::FLAMELET: + AddVolumeOutput("PROGVAR", "Progress_Variable", "SOLUTION", "Progress variable"); + AddVolumeOutput("ENTHALPY", "Total_Enthalpy", "SOLUTION", "Total enthalpy"); + /*--- auxiliary species ---*/ + for (unsigned short iReactant=0; iReactantGetNUserScalars(); iReactant++) { + string species_name = config->GetUserScalarName(iReactant); + AddVolumeOutput(species_name, species_name, "SOLUTION", species_name + "Mass fraction solution"); + } + + AddVolumeOutput("SOURCE_PROGVAR", "Source_Progress_Variable", "SOURCE", "Source Progress Variable"); + /*--- no source term for enthalpy ---*/ + /*--- auxiliary species source terms ---*/ + for (unsigned short iReactant=0; iReactantGetNUserScalars(); iReactant++) { + string species_name = config->GetUserScalarName(iReactant); + AddVolumeOutput("SOURCE_" + species_name, "Source_" + species_name, "SOURCE", "Source " + species_name); + } + + for (int i_lookup = 0; i_lookup < config->GetNLookups(); ++i_lookup) { + if (config->GetLUTLookupName(i_lookup) != "NULL") { + string strname1 = "lookup_" + config->GetLUTLookupName(i_lookup); + AddVolumeOutput(config->GetLUTLookupName(i_lookup), strname1,"LOOKUP", config->GetLUTLookupName(i_lookup)); + } + } + + AddVolumeOutput("TABLE_MISSES" , "Table_misses" , "SOLUTION", "Lookup table misses"); + break; + case SPECIES_MODEL::NONE: + break; } } @@ -1131,6 +1274,24 @@ void CFlowOutput::SetVolumeOutputFields_ScalarResidual(const CConfig* config) { break; } + switch (config->GetKind_Species_Model()) { + case SPECIES_MODEL::SPECIES_TRANSPORT: + for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) + AddVolumeOutput("RES_SPECIES_" + std::to_string(iVar), "Residual_Species_" + std::to_string(iVar), "RESIDUAL", "Residual of the transported species " + std::to_string(iVar)); + break; + case SPECIES_MODEL::FLAMELET: + AddVolumeOutput("RES_PROGVAR", "Residual_Progress_Variable", "SOLUTION", "Residual of progress variable"); + AddVolumeOutput("RES_ENTHALPY", "Residual_Total_Enthalpy", "SOLUTION", "Residual of total enthalpy"); + /*--- residuals for auxiliary species transport equations ---*/ + for (unsigned short iReactant=0; iReactantGetNUserScalars(); iReactant++){ + string species_name = config->GetUserScalarName(iReactant); + AddVolumeOutput("RES_" + species_name, "Residual_" + species_name, "RESIDUAL", "Residual of the " + species_name + " equation"); + } + break; + case SPECIES_MODEL::NONE: + break; + } + switch (config->GetKind_Trans_Model()) { case TURB_TRANS_MODEL::LM: AddVolumeOutput("RES_INTERMITTENCY", "Residual_LM_intermittency", "RESIDUAL", "Residual of LM intermittency"); @@ -1140,11 +1301,6 @@ void CFlowOutput::SetVolumeOutputFields_ScalarResidual(const CConfig* config) { case TURB_TRANS_MODEL::NONE: break; } - - if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { - for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) - AddVolumeOutput("RES_SPECIES_" + std::to_string(iVar), "Residual_Species_" + std::to_string(iVar), "RESIDUAL", "Residual of the transported species " + std::to_string(iVar)); - } } void CFlowOutput::SetVolumeOutputFields_ScalarLimiter(const CConfig* config) { @@ -1164,10 +1320,23 @@ void CFlowOutput::SetVolumeOutputFields_ScalarLimiter(const CConfig* config) { } } - if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { - if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) { - for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) - AddVolumeOutput("LIMITER_SPECIES_" + std::to_string(iVar), "Limiter_Species_" + std::to_string(iVar), "LIMITER", "Limiter value of the transported species " + std::to_string(iVar)); + if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) { + switch (config->GetKind_Species_Model()) { + case SPECIES_MODEL::SPECIES_TRANSPORT: + for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) + AddVolumeOutput("LIMITER_SPECIES_" + std::to_string(iVar), "Limiter_Species_" + std::to_string(iVar), "LIMITER", "Limiter value of the transported species " + std::to_string(iVar)); + break; + case SPECIES_MODEL::FLAMELET: + AddVolumeOutput("LIMITER_PROGVAR", "Limiter_Progress_Variable", "SOLUTION", "Limiter of progress variable"); + AddVolumeOutput("LIMITER_ENTHALPY", "Limiter_Total_Enthalpy", "SOLUTION", "Limiter of total enthalpy"); + /*--- limiter for auxiliary species transport ---*/ + for (unsigned short iReactant=0; iReactantGetNUserScalars(); iReactant++) { + string species_name = config->GetUserScalarName(iReactant); + AddVolumeOutput("LIMITER_" + species_name, "LIMITER_" + species_name, "LIMITER", "Limiter value for the " + species_name + " equation"); + } + break; + case SPECIES_MODEL::NONE: + break; } } @@ -1270,15 +1439,55 @@ void CFlowOutput::LoadVolumeData_Scalar(const CConfig* config, const CSolver* co SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint)); } - if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { - const auto Node_Species = solver[SPECIES_SOL]->GetNodes(); + switch (config->GetKind_Species_Model()) { + case SPECIES_MODEL::SPECIES_TRANSPORT: { + const auto Node_Species = solver[SPECIES_SOL]->GetNodes(); + for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { + SetVolumeOutputValue("SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetSolution(iPoint, iVar)); + SetVolumeOutputValue("RES_SPECIES_" + std::to_string(iVar), iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, iVar)); + if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) + SetVolumeOutputValue("LIMITER_SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetLimiter(iPoint, iVar)); + } + break; + } - for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - SetVolumeOutputValue("SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetSolution(iPoint, iVar)); - SetVolumeOutputValue("RES_SPECIES_" + std::to_string(iVar), iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, iVar)); - if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) - SetVolumeOutputValue("LIMITER_SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetLimiter(iPoint, iVar)); + case SPECIES_MODEL::FLAMELET: { + const auto Node_Species = solver[SPECIES_SOL]->GetNodes(); + + SetVolumeOutputValue("PROGVAR", iPoint, Node_Species->GetSolution(iPoint, I_PROGVAR)); + SetVolumeOutputValue("ENTHALPY", iPoint, Node_Species->GetSolution(iPoint, I_ENTH)); + SetVolumeOutputValue("SOURCE_PROGVAR", iPoint, Node_Species->GetScalarSources(iPoint, I_PROGVAR)); + SetVolumeOutputValue("RES_PROGVAR", iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, I_PROGVAR)); + SetVolumeOutputValue("RES_ENTHALPY", iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, I_ENTH)); + SetVolumeOutputValue("TABLE_MISSES" , iPoint, (su2double)Node_Species->GetInsideTable(iPoint)); + + /*--- auxiliary species transport equations ---*/ + for (unsigned short i_scalar=0; i_scalarGetNUserScalars(); i_scalar++) { + string scalar_name = config->GetUserScalarName(i_scalar); + SetVolumeOutputValue(scalar_name, iPoint, Node_Species->GetSolution(iPoint, config->GetNControlVars() + i_scalar)); + SetVolumeOutputValue("SOURCE_" + scalar_name, iPoint, Node_Species->GetScalarSources(iPoint, config->GetNControlVars() + i_scalar)); + SetVolumeOutputValue("RES_" + scalar_name, iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, config->GetNControlVars() + i_scalar)); + } + + if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) { + SetVolumeOutputValue("LIMITER_PROGVAR", iPoint, Node_Species->GetLimiter(iPoint, I_PROGVAR)); + SetVolumeOutputValue("LIMITER_ENTHALPY", iPoint, Node_Species->GetLimiter(iPoint, I_ENTH)); + /*--- limiter for auxiliary species transport equations ---*/ + for (unsigned short i_scalar=0; i_scalarGetNUserScalars(); i_scalar++) { + string scalar_name = config->GetUserScalarName(i_scalar); + SetVolumeOutputValue("LIMITER_" + scalar_name, iPoint, Node_Species->GetLimiter(iPoint, config->GetNControlVars() + i_scalar)); + } + } + + /*--- variables that we look up from the LUT ---*/ + for (int i_lookup = 0; i_lookup < config->GetNLookups(); ++i_lookup) { + if (config->GetLUTLookupName(i_lookup)!="NULL") + SetVolumeOutputValue(config->GetLUTLookupName(i_lookup), iPoint, Node_Species->GetScalarLookups(iPoint, i_lookup)); + } + + break; } + case SPECIES_MODEL::NONE: break; } } @@ -2428,6 +2637,10 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo << config->GetTemperature_Critical() / config->GetTemperature_Ref() << "\n"; break; + case FLUID_FLAMELET: + file << "Fluid Model: FLAMELET \n"; + break; + case COOLPROP: { CCoolProp auxFluidModel(config->GetFluid_Name()); file << "Fluid Model: CoolProp library \n"; @@ -2747,6 +2960,12 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo else file << " psf.\n"; break; + case FLUID_FLAMELET: + file << "Fluid model: FLUID_FLAMELET \n"; + if (si_units) file << " Pa.\n"; + else file << " psf.\n"; + break; + case INC_IDEAL_GAS_POLY: file << "Fluid Model: INC_IDEAL_GAS_POLY \n"; file << "Variable density incompressible flow using ideal gas law.\n"; @@ -2781,6 +3000,13 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo file << "Laminar Viscosity (non-dim): " << config->GetMu_ConstantND() << "\n"; break; + case VISCOSITYMODEL::FLAMELET: + file << "Viscosity Model: FLAMELET \n"; + if (si_units) file << " N.s/m^2.\n"; + else file << " lbf.s/ft^2.\n"; + file << "Laminar Viscosity (non-dim): " << config->GetMu_ConstantND() << "\n"; + break; + case VISCOSITYMODEL::COOLPROP: file << "Viscosity Model: CoolProp \n"; break; @@ -2834,6 +3060,12 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo file << "Conductivity Model: COOLPROP \n"; break; + case CONDUCTIVITYMODEL::FLAMELET: + file << "Conductivity Model: FLAMELET \n"; + file << "Molecular Conductivity units: " << " W/m^2.K.\n"; + file << "Molecular Conductivity (non-dim): " << config->GetThermal_Conductivity_ConstantND() << "\n"; + break; + case CONDUCTIVITYMODEL::POLYNOMIAL: file << "Viscosity Model: POLYNOMIAL \n"; file << "Kt(T) polynomial coefficients: \n ("; diff --git a/SU2_CFD/src/output/output_structure_legacy.cpp b/SU2_CFD/src/output/output_structure_legacy.cpp index cdf6a468c74b..11a4be23942a 100644 --- a/SU2_CFD/src/output/output_structure_legacy.cpp +++ b/SU2_CFD/src/output/output_structure_legacy.cpp @@ -3319,6 +3319,13 @@ void COutputLegacy::SpecialOutput_ForcesBreakdown(CSolver *****solver, CGeometry Breakdown_file << "Laminar Viscosity (non-dim): " << config[val_iZone]->GetMu_ConstantND()<< "\n"; break; + case VISCOSITYMODEL::FLAMELET: + Breakdown_file << "Viscosity Model: FLAMELET "<< "\n"; + if (config[val_iZone]->GetSystemMeasurements() == SI) Breakdown_file << " N.s/m^2." << "\n"; + else if (config[val_iZone]->GetSystemMeasurements() == US) Breakdown_file << " lbf.s/ft^2." << "\n"; + Breakdown_file << "Laminar Viscosity (non-dim): " << config[val_iZone]->GetMu_ConstantND()<< "\n"; + break; + case VISCOSITYMODEL::SUTHERLAND: Breakdown_file << "Viscosity Model: SUTHERLAND "<< "\n"; Breakdown_file << "Ref. Laminar Viscosity: " << config[val_iZone]->GetMu_Ref(); @@ -3371,6 +3378,12 @@ void COutputLegacy::SpecialOutput_ForcesBreakdown(CSolver *****solver, CGeometry Breakdown_file << "Molecular Conductivity (non-dim): " << config[val_iZone]->GetThermal_Conductivity_ConstantND()<< "\n"; break; + case CONDUCTIVITYMODEL::FLAMELET: + Breakdown_file << "Conductivity Model: FLAMELET "<< "\n"; + Breakdown_file << "Molecular Conductivity units: " << " W/m^2.K." << "\n"; + Breakdown_file << "Molecular Conductivity (non-dim): " << config[val_iZone]->GetThermal_Conductivity_ConstantND()<< "\n"; + break; + case CONDUCTIVITYMODEL::COOLPROP: Breakdown_file << "Conductivity Model: COOLPROP "<< "\n"; break; diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 32672e87eef3..adbdc9099dc7 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -24,7 +24,7 @@ * You should have received a copy of the GNU Lesser General Public * License along with SU2. If not, see . */ - + #define ENABLE_MAPS #include "../include/drivers/CDriver.hpp" #include "../include/drivers/CSinglezoneDriver.hpp" @@ -35,7 +35,7 @@ void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geom int rank = MASTER_NODE; SU2_MPI::Comm_rank(SU2_MPI::GetComm(), &rank); - /* --- Initialize boundary conditions customization, this is achieve through the Python wrapper --- */ + /* --- Initialize boundary conditions customization, this is achieved through the Python wrapper --- */ for(iZone=0; iZone < nZone; iZone++){ if (config[iZone]->GetnMarker_PyCustom() > 0){ @@ -55,6 +55,11 @@ void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geom solver[iZone][INST_0][MESH_0][FLOW_SOL]->UpdateCustomBoundaryConditions(geometry[iZone][INST_0], config[iZone]); } } + + if (config[iZone]->GetInitial_PyCustom() > 0){ + if (rank == MASTER_NODE) cout << endl << "----------------- Python Initialization Preprocessing ( Zone "<< iZone <<" ) -----------------" << endl; + } + } } @@ -176,6 +181,24 @@ unsigned long CDriver::GetNumberVertices(unsigned short iMarker) const { } +/* coordinates of the point on the mesh*/ +vector CDriver::GetCoords(unsigned long iPoint, unsigned short iMesh) const { + vector coord(3,0.0); + vector coord_passive(3, 0.0); + + for (auto iDim = 0 ; iDim < nDim ; iDim++){ + coord[iDim] = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint,iDim); + } + + coord_passive[0] = SU2_TYPE::GetValue(coord[0]); + coord_passive[1] = SU2_TYPE::GetValue(coord[1]); + coord_passive[2] = SU2_TYPE::GetValue(coord[2]); + + return coord_passive; +} + +unsigned long CDriver::GetNumberVertices() const { return geometry_container[ZONE_0][INST_0][MESH_0]->GetnPoint(); } + unsigned long CDriver::GetNumberHaloVertices(unsigned short iMarker) const { unsigned long nHaloVertices, iVertex, iPoint; @@ -211,6 +234,11 @@ bool CDriver::IsAHaloNode(unsigned short iMarker, unsigned long iVertex) const { } +/*--- Get the index to the solver ---*/ +unsigned long CDriver::GetSolverIndex(string solverName) const { + return static_cast(SolverType_Map.find(solverName)->second); +} + vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) const { vector coord(3,0.0); @@ -268,6 +296,11 @@ unsigned long CDriver::GetnTimeIter() const { return config_container[ZONE_0]->GetnTime_Iter(); } +// number of points on the mesh +unsigned long CDriver::GetnPoints(unsigned short iMesh) const { + return geometry_container[ZONE_0][INST_0][iMesh]->GetnPointDomain(); +} + unsigned long CDriver::GetTime_Iter() const{ return TimeIter; @@ -556,6 +589,66 @@ void CDriver::SetHeatSource_Position(passivedouble alpha, passivedouble pos_x, p } +void CDriver::SetStates(vector> values) { + const auto nPoint = GetNumberVertices(); + + if (values.size() != nPoint) { + SU2_MPI::Error("Invalid number of vertices!", CURRENT_FUNCTION); + } + + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + SetStates(iPoint, values[iPoint]); + } +} + +void CDriver::SetStates(unsigned long iPoint, vector values) { + const auto nVar = GetNumberStateVariables(); + + if (values.size() != nVar) { + SU2_MPI::Error("Invalid number of variables!", CURRENT_FUNCTION); + } + + for (auto iVar = 0u; iVar < nVar; iVar++) { + solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->SetSolution(iPoint, iVar, values[iVar]); + } +} + +void CDriver::SetStates(vector> values, const int SOLVER) { + const auto nPoint = GetNumberVertices(); + + if (values.size() != nPoint) { + SU2_MPI::Error("Invalid number of vertices!", CURRENT_FUNCTION); + } + + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + SetStates(iPoint, values[iPoint]); + } +} + +void CDriver::SetStates(unsigned long iPoint, vector values, const int SOLVER) { + const auto nVar = GetNumberStateVariables(SOLVER); + if (values.size() != nVar) { + SU2_MPI::Error("Invalid number of variables!", CURRENT_FUNCTION); + } + + for (auto iVar = 0u; iVar < nVar; iVar++) { + solver_container[ZONE_0][INST_0][MESH_0][SOLVER]->GetNodes()->SetSolution(iPoint, iVar, values[iVar]); + } +} + +unsigned long CDriver::GetNumberStateVariables() const { + if (!config_container[ZONE_0]->GetFluidProblem()) { + SU2_MPI::Error("Flow solver is not defined!", CURRENT_FUNCTION); + } + + return solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetnVar(); +} + +unsigned long CDriver::GetNumberStateVariables(const int SOLVER) const { + + return solver_container[ZONE_0][INST_0][MESH_0][SOLVER]->GetnVar(); +} + void CDriver::SetInlet_Angle(unsigned short iMarker, passivedouble alpha){ su2double alpha_rad = alpha * PI_NUMBER/180.0; @@ -893,3 +986,6 @@ vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long return FlowLoad_passive; } + + #undef ENABLE_MAPS + diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index f44f76817fa1..c1b31b9726e2 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1093,6 +1093,7 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes if (viscous) { GetFluidModel()->SetLaminarViscosityModel(config); GetFluidModel()->SetThermalConductivityModel(config); + GetFluidModel()->SetMassDiffusivityModel(config); } } diff --git a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp index 7acc893e9861..18b5d6f8b517 100644 --- a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp +++ b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp @@ -1090,6 +1090,7 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig *config, if (viscous) { FluidModel->SetLaminarViscosityModel(config); FluidModel->SetThermalConductivityModel(config); + FluidModel->SetMassDiffusivityModel(config); // nijso: TODO, needs to be tested } if (tkeNeeded) { Energy_FreeStreamND += Tke_FreeStreamND; }; config->SetEnergy_FreeStreamND(Energy_FreeStreamND); diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 1749ae096f22..7f7fb64c796a 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -34,6 +34,7 @@ #include "../../include/limiters/CLimiterDetails.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../include/fluid/CFluidScalar.hpp" +#include "../../include/fluid/CFluidFlamelet.hpp" #include "../../include/fluid/CFluidModel.hpp" @@ -262,6 +263,8 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i bool energy = config->GetEnergy_Equation(); bool boussinesq = (config->GetKind_DensityModel() == INC_DENSITYMODEL::BOUSSINESQ); + bool redefine_fluid_model = true; + /*--- Compute dimensional free-stream values. ---*/ Density_FreeStream = config->GetInc_Density_Init(); config->SetDensity_FreeStream(Density_FreeStream); @@ -328,6 +331,18 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i auxFluidModel->SetTDState_T(Temperature_FreeStream, config->GetSpecies_Init()); break; + case FLUID_FLAMELET: + + config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT / (config->GetMolecular_Weight() / 1000.0)); + Pressure_Thermodynamic = Density_FreeStream * Temperature_FreeStream * config->GetGas_Constant(); + auxFluidModel = new CFluidFlamelet(config, Pressure_Thermodynamic); + auxFluidModel->SetTDState_T(Temperature_FreeStream, config->GetSpecies_Init()); + config->SetPressure_Thermodynamic(Pressure_Thermodynamic); + + /*--- flamelet fluid model does not need to be redefined after nondimensionalization ---*/ + redefine_fluid_model = false; + break; + default: SU2_MPI::Error("Fluid model not implemented for incompressible solver.", CURRENT_FUNCTION); @@ -452,7 +467,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i /*--- Delete the original (dimensional) FluidModel object. No fluid is used for inscompressible cases. ---*/ - delete auxFluidModel; + if(redefine_fluid_model) delete auxFluidModel; /*--- Create one final fluid model object per OpenMP thread to be able to use them in parallel. * GetFluidModel() should be used to automatically access the "right" object of each thread. ---*/ @@ -478,6 +493,11 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i fluidModel->SetTDState_T(Temperature_FreeStreamND, config->GetSpecies_Init()); break; + case FLUID_FLAMELET: + fluidModel = auxFluidModel; + fluidModel->SetTDState_T(Temperature_FreeStreamND, config->GetSpecies_Init()); + break; + case INC_IDEAL_GAS_POLY: fluidModel = new CIncIdealGasPolynomial(Gas_ConstantND, Pressure_ThermodynamicND); if (viscous) { @@ -510,7 +530,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i fluidModel->SetLaminarViscosityModel(config); fluidModel->SetThermalConductivityModel(config); - + fluidModel->SetMassDiffusivityModel(config); } } @@ -625,6 +645,15 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i NonDimTable.PrintFooter(); break; + case VISCOSITYMODEL::FLAMELET: + ModelTable << "FLAMELET"; + if (config->GetSystemMeasurements() == SI) Unit << "N.s/m^2"; + else if (config->GetSystemMeasurements() == US) Unit << "lbf.s/ft^2"; + NonDimTable << "Viscosity" << "--" << "--" << Unit.str() << config->GetMu_ConstantND(); + Unit.str(""); + NonDimTable.PrintFooter(); + break; + case VISCOSITYMODEL::COOLPROP: ModelTable << "COOLPROP_VISCOSITY"; if (config->GetSystemMeasurements() == SI) Unit << "N.s/m^2"; @@ -682,6 +711,13 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i NonDimTable.PrintFooter(); break; + case CONDUCTIVITYMODEL::FLAMELET: + ModelTable << "FLAMELET"; + Unit << "W/m^2.K"; + NonDimTable << "Molecular Cond." << "--" << "--" << Unit.str() << config->GetThermal_Conductivity_ConstantND(); + Unit.str(""); + break; + case CONDUCTIVITYMODEL::COOLPROP: ModelTable << "COOLPROP"; Unit << "W/m^2.K"; @@ -759,6 +795,23 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i NonDimTable.PrintFooter(); break; + case FLUID_FLAMELET: + ModelTable << "FLAMELET"; + Unit << "N.m/kg.K"; + NonDimTable << "Spec. Heat (Cp)" << "--" << "--" << Unit.str() << config->GetSpecific_Heat_CpND(); + Unit.str(""); + Unit << "g/mol"; + NonDimTable << "Molecular weight" << "--" << "--" << Unit.str() << config->GetMolecular_Weight(); + Unit.str(""); + Unit << "N.m/kg.K"; + NonDimTable << "Gas Constant" << "--" << config->GetGas_Constant_Ref() << Unit.str() << config->GetGas_ConstantND(); + Unit.str(""); + Unit << "Pa"; + NonDimTable << "Therm. Pressure" << config->GetPressure_Thermodynamic() << config->GetPressure_Ref() << Unit.str() << config->GetPressure_ThermodynamicND(); + Unit.str(""); + NonDimTable.PrintFooter(); + break; + case INC_IDEAL_GAS_POLY: ModelTable << "INC_IDEAL_GAS_POLY"; Unit.str(""); @@ -2904,8 +2957,6 @@ void CIncEulerSolver::GetOutlet_Properties(CGeometry *geometry, CConfig *config, if (geometry->nodes->GetDomain(iPoint)) { - V_outlet = nodes->GetPrimitive(iPoint); - geometry->vertex[iMarker][iVertex]->GetNormal(Vector); su2double AxiFactor = 1.0; @@ -2921,7 +2972,9 @@ void CIncEulerSolver::GetOutlet_Properties(CGeometry *geometry, CConfig *config, } } - Density = V_outlet[prim_idx.Density()]; + V_outlet = nodes->GetPrimitive(iPoint); + + Density = V_outlet[prim_idx.Density()]; Velocity2 = 0.0; Area = 0.0; MassFlow = 0.0; @@ -2931,6 +2984,7 @@ void CIncEulerSolver::GetOutlet_Properties(CGeometry *geometry, CConfig *config, Velocity2 += Velocity[iDim] * Velocity[iDim]; MassFlow += Vector[iDim] * AxiFactor * Density * Velocity[iDim]; } + Area = sqrt (Area); Outlet_MassFlow[iMarker] += MassFlow; diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 357b43dee867..85cabf300ece 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -3659,6 +3659,14 @@ void CSolver::LoadInletProfile(CGeometry **geometry, columnValue << config->GetInlet_SpeciesVal(Marker_Tag)[iVar] << "\t"; } break; + case SPECIES_MODEL::FLAMELET: + /*--- 2-equation flamelet model ---*/ + columnName << "PROGVAR " << setw(24) << "ENTHALPY " << setw(24); + columnValue << config->GetInlet_SpeciesVal(Marker_Tag)[0] << "\t" << config->GetInlet_SpeciesVal(Marker_Tag)[1]<<"\t"; + /*--- auxiliary species transport equations ---*/ + columnName << "CO " << setw(24) << "NOX " << setw(24); + columnValue << config->GetInlet_SpeciesVal(Marker_Tag)[2] << "\t" << config->GetInlet_SpeciesVal(Marker_Tag)[3]<<"\t"; + break; } columnNames.push_back(columnName.str()); diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index 02843124a4f2..5ea071ee99da 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -52,6 +52,7 @@ #include "../../include/solvers/CBaselineSolver_FEM.hpp" #include "../../include/solvers/CRadP1Solver.hpp" #include "../../include/solvers/CSpeciesSolver.hpp" +#include "../../include/solvers/CSpeciesFlameletSolver.hpp" map CSolverFactory::allocatedSolvers; @@ -394,13 +395,25 @@ CSolver* CSolverFactory::CreateSpeciesSolver(CSolver **solver, CGeometry *geomet CSolver *speciesSolver = nullptr; - if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { - if (adjoint){ - speciesSolver = new CDiscAdjSolver(geometry, config, solver[SPECIES_SOL], RUNTIME_SPECIES_SYS, iMGLevel); - } else { - speciesSolver = new CSpeciesSolver(geometry, config, iMGLevel); - } + switch (config->GetKind_Species_Model()) { + case SPECIES_MODEL::SPECIES_TRANSPORT: + if (adjoint) { + speciesSolver = new CDiscAdjSolver(geometry, config, solver[SPECIES_SOL], RUNTIME_SPECIES_SYS, iMGLevel); + } else { + speciesSolver = new CSpeciesSolver(geometry, config, iMGLevel); + } + break; + case SPECIES_MODEL::FLAMELET: + if (adjoint) { + speciesSolver = new CDiscAdjSolver(geometry, config, solver[SPECIES_SOL], RUNTIME_SPECIES_SYS, iMGLevel); + } else { + speciesSolver = new CSpeciesFlameletSolver(geometry, config, iMGLevel); + } + break; + case SPECIES_MODEL::NONE: + break; } + return speciesSolver; } @@ -411,9 +424,9 @@ CSolver* CSolverFactory::CreateHeatSolver(CSolver **solver, CGeometry *geometry, /*--- Only allocate a heat solver if it should run standalone * or if the weakly coupled heat solver is enabled and no energy equation is included ---*/ - if ((config->GetWeakly_Coupled_Heat() && !config->GetEnergy_Equation()) || config->GetHeatProblem()){ - if (adjoint){ - if (config->GetDiscrete_Adjoint()){ + if ((config->GetWeakly_Coupled_Heat() && !config->GetEnergy_Equation()) || config->GetHeatProblem()) { + if (adjoint) { + if (config->GetDiscrete_Adjoint()) { heatSolver = new CDiscAdjSolver(geometry, config, solver[HEAT_SOL], RUNTIME_HEAT_SYS, iMGLevel); } else { @@ -424,24 +437,24 @@ CSolver* CSolverFactory::CreateHeatSolver(CSolver **solver, CGeometry *geometry, heatSolver = new CHeatSolver(geometry, config, iMGLevel); } } - return heatSolver; + return heatSolver; } CSolver* CSolverFactory::CreateMeshSolver(CSolver **solver, CGeometry *geometry, CConfig *config, int iMGLevel, bool adjoint){ CSolver *meshSolver = nullptr; - if (config->GetDeform_Mesh() && iMGLevel == MESH_0){ - if (!adjoint){ + if (config->GetDeform_Mesh() && iMGLevel == MESH_0) { + if (!adjoint) { meshSolver = new CMeshSolver(geometry, config); } - if (adjoint && config->GetDiscrete_Adjoint()){ + if (adjoint && config->GetDiscrete_Adjoint()) { meshSolver = new CDiscAdjMeshSolver(geometry, config, solver[MESH_SOL]); } } - return meshSolver; + return meshSolver; } CSolver* CSolverFactory::CreateDGSolver(SUB_SOLVER_TYPE kindDGSolver, CGeometry *geometry, CConfig *config, int iMGLevel){ diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp new file mode 100644 index 000000000000..98f38fcab4f0 --- /dev/null +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -0,0 +1,817 @@ +/*! + * \file CSpeciesFlameletSolver.cpp + * \brief Main subroutines of CSpeciesFlameletSolver class + * \author D. Mayer, T. Economon, N. Beishuizen + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/solvers/CSpeciesFlameletSolver.hpp" + +#include "../../../Common/include/parallelization/omp_structure.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../include/fluid/CFluidFlamelet.hpp" +#include "../../include/solvers/CSpeciesSolver.hpp" +#include "../../include/variables/CFlowVariable.hpp" +#include "../../include/variables/CSpeciesFlameletVariable.hpp" + +CSpeciesFlameletSolver::CSpeciesFlameletSolver(CGeometry* geometry, CConfig* config, unsigned short iMesh) + : CSpeciesSolver(geometry, config, true) { + /*--- Store if an implicit scheme is used, for use during periodic boundary conditions. ---*/ + SetImplicitPeriodic(config->GetKind_TimeIntScheme_Species() == EULER_IMPLICIT); + + /*--- Dimension of the problem. ---*/ + nVar = config->GetNScalars(); + nPrimVar = nVar; + + if (nVar > MAXNVAR) + SU2_MPI::Error("Increase static array size MAXNVAR for CSpeciesVariable and proceed.", CURRENT_FUNCTION); + + nPoint = geometry->GetnPoint(); + nPointDomain = geometry->GetnPointDomain(); + + /*--- Initialize nVarGrad for deallocation ---*/ + + nVarGrad = nVar; + + /*--- Define geometry constants in the solver structure ---*/ + + nDim = geometry->GetnDim(); + + /*--- Define some auxiliary vector related with the solution ---*/ + Solution = new su2double[nVar]; + Solution_i = new su2double[nVar]; + Solution_j = new su2double[nVar]; + + /*--- Allocates a 3D array with variable "middle" sizes and init to 0. ---*/ + + auto Alloc3D = [](unsigned long M, const vector& N, unsigned long P, vector& X) { + X.resize(M); + for (unsigned long i = 0; i < M; ++i) X[i].resize(N[i], P) = su2double(0.0); + }; + + /*--- Store the values of the temperature and the heat flux density at the boundaries, + used for coupling with a solid donor cell ---*/ + constexpr auto n_conjugate_var = 4u; + + Alloc3D(nMarker, nVertex, n_conjugate_var, conjugate_var); + for (auto& x : conjugate_var) x = config->GetTemperature_FreeStreamND(); + + /*--- Single grid simulation ---*/ + + if (iMesh == MESH_0 || config->GetMGCycle() == FULLMG_CYCLE) { + /*--- Define some auxiliary vector related with the residual ---*/ + + Residual_RMS.resize(nVar, 0.0); + Residual_Max.resize(nVar, 0.0); + Point_Max.resize(nVar, 0); + Point_Max_Coord.resize(nVar, nDim) = su2double(0.0); + + /*--- Initialize the BGS residuals in multizone problems. ---*/ + if (config->GetMultizone_Problem()) { + Residual_BGS.resize(nVar, 0.0); + Residual_Max_BGS.resize(nVar, 0.0); + Point_Max_BGS.resize(nVar, 0); + Point_Max_Coord_BGS.resize(nVar, nDim) = su2double(0.0); + } + + /*--- Initialization of the structure of the whole Jacobian ---*/ + + if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (flamelet model)." << endl; + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); + LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); + LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); + System.SetxIsZero(true); + + if (ReducerStrategy) EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + } + + /*--- Initialize lower and upper limits---*/ + + if (config->GetSpecies_Clipping()) { + for (auto iVar = 0u; iVar < nVar; iVar++) { + lowerlimit[iVar] = config->GetSpecies_Clipping_Min(iVar); + upperlimit[iVar] = config->GetSpecies_Clipping_Max(iVar); + } + } else { + for (auto iVar = 0u; iVar < nVar; iVar++) { + /*--- we fix the lower limit to 0 ---*/ + lowerlimit[iVar] = -1.0e15; + upperlimit[iVar] = 1.0e15; + } + } + + /*--- Scalar variable state at the far-field. ---*/ + + for (auto iVar = 0u; iVar < nVar; iVar++) { + Solution_Inf[iVar] = config->GetSpecies_Init()[iVar]; + } + + /*--- Initialize the solution to the far-field state everywhere. ---*/ + + nodes = new CSpeciesFlameletVariable(Solution_Inf, nPoint, nDim, nVar, config); + SetBaseClassPointerToNodes(); + + /*--- MPI solution ---*/ + + InitiateComms(geometry, config, SOLUTION); + CompleteComms(geometry, config, SOLUTION); + + /*--- Set the column number for species in inlet-files. + * e.g. Coords(nDim), Temp(1), VelMag(1), Normal(nDim), Turb(1 or 2), Species(arbitrary) ---*/ + Inlet_Position = nDim + 2 + nDim + config->GetnTurbVar(); + + /*-- Allocation of inlet-values. Will be filled either by an inlet files, + * or uniformly by a uniform boundary condition. ---*/ + + Inlet_SpeciesVars.resize(nMarker); + for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) { + Inlet_SpeciesVars[iMarker].resize(nVertex[iMarker], nVar); + for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; ++iVertex) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Inlet_SpeciesVars[iMarker](iVertex, iVar) = Solution_Inf[iVar]; + } + } + } + + /*--- Store the initial CFL number for all grid points. ---*/ + + const su2double CFL = config->GetCFL(MGLevel) * config->GetCFLRedCoeff_Species(); + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0u; iPoint < nPoint; iPoint++) { + nodes->SetLocalCFL(iPoint, CFL); + } + END_SU2_OMP_FOR + Min_CFL_Local = CFL; + Max_CFL_Local = CFL; + Avg_CFL_Local = CFL; + + /*--- Add the solver name (max 8 characters) ---*/ + SolverName = "FLAMELET"; +} + +void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, + unsigned short iMesh, unsigned short iRKStep, + unsigned short RunTime_EqSystem, bool Output) { + unsigned long n_not_in_domain = 0; + unsigned long global_table_misses = 0; + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes(); + + for (auto i_point = 0u; i_point < nPoint; i_point++) { + CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); + su2double* scalars = nodes->GetSolution(i_point); + + /*--- Compute scalar source terms ---*/ + unsigned long exit_code = fluid_model_local->SetScalarSources(scalars); + unsigned short inside = exit_code; + nodes->SetInsideTable(i_point, inside); + n_not_in_domain += exit_code; + + /*--- Get lookup scalars ---*/ + fluid_model_local->SetScalarLookups(scalars); + for(auto i_lookup=0u; i_lookupGetNLookups(); i_lookup++){ + nodes->SetLookupScalar(i_point, fluid_model_local->GetScalarLookups(i_lookup), i_lookup); + } + + for(auto i_scalar=0u; i_scalar < nVar; i_scalar++) + nodes->SetScalarSource(i_point, i_scalar, fluid_model_local->GetScalarSources(i_scalar)); + + su2double 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); + } + + if (!Output) LinSysRes.SetBlock_Zero(i_point); + } + + if (config->GetComm_Level() == COMM_FULL) { + SU2_MPI::Reduce(&n_not_in_domain, &global_table_misses, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, + SU2_MPI::GetComm()); + if (rank == MASTER_NODE) { + SetNTableMisses(global_table_misses); + } + } + + /*--- Clear residual and system matrix, not needed for + * reducer strategy as we write over the entire matrix. ---*/ + if (!ReducerStrategy && !Output) { + LinSysRes.SetValZero(); + if (implicit) + Jacobian.SetValZero(); + else { + SU2_OMP_BARRIER + } + } + + /*--- Upwind second order reconstruction and gradients ---*/ + + if (config->GetReconstructionGradientRequired()) { + if (config->GetKind_Gradient_Method_Recon() == GREEN_GAUSS) SetSolution_Gradient_GG(geometry, config, true); + if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) SetSolution_Gradient_LS(geometry, config, true); + if (config->GetKind_Gradient_Method_Recon() == WEIGHTED_LEAST_SQUARES) + SetSolution_Gradient_LS(geometry, config, true); + } + + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) SetSolution_Gradient_GG(geometry, config); + + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) SetSolution_Gradient_LS(geometry, config); + +} + +void CSpeciesFlameletSolver::Postprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, + unsigned short iMesh) { + /*--- your postprocessing goes here ---*/ +} + +void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver*** solver_container, CConfig* config, + unsigned long ExtIter) { + bool Restart = (config->GetRestart() || config->GetRestart_Flow()); + + /*--- do not use initial condition when custom python is active ---*/ + if (config->GetInitial_PyCustom()) { + if (rank == MASTER_NODE) cout << "Initialization through custom python function." << endl; + return; + } + + if ((!Restart) && ExtIter == 0) { + if (rank == MASTER_NODE) { + cout << "Initializing progress variable and total enthalpy (using temperature)" << endl; + } + + su2double* scalar_init = new su2double[nVar]; + su2double* flame_offset = config->GetFlameOffset(); + su2double* flame_normal = config->GetFlameNormal(); + + su2double prog_burnt; + su2double prog_unburnt = 0.0; + su2double flame_thickness = config->GetFlameThickness(); + su2double burnt_thickness = config->GetFlameBurntThickness(); + su2double flamenorm = GeometryToolbox::Norm(nDim, flame_normal); + + su2double temp_inlet = config->GetInc_Temperature_Init(); + su2double prog_inlet = config->GetSpecies_Init()[I_PROGVAR]; + su2double enth_inlet = config->GetSpecies_Init()[I_ENTH]; + if (rank == MASTER_NODE) { + cout << "initial condition: T = " << temp_inlet << endl; + cout << "initial condition: c = " << prog_inlet << endl; + cout << "initial condition: h = " << enth_inlet << endl; + } + + su2double point_loc; + + CFluidModel* fluid_model_local; + + vector look_up_tags; + vector look_up_data; + string name_enth = config->GetLUTScalarName(I_ENTH); + string name_prog = config->GetLUTScalarName(I_PROGVAR); + + unsigned long n_not_iterated_local = 0; + unsigned long n_not_in_domain_local = 0; + unsigned long n_points_unburnt_local = 0; + unsigned long n_points_burnt_local = 0; + unsigned long n_points_flame_local = 0; + unsigned long n_not_iterated_global; + unsigned long n_not_in_domain_global; + unsigned long n_points_burnt_global; + unsigned long n_points_flame_global; + unsigned long n_points_unburnt_global; + + for (unsigned long i_mesh = 0; i_mesh <= config->GetnMGLevels(); i_mesh++) { + fluid_model_local = solver_container[i_mesh][FLOW_SOL]->GetFluidModel(); + + prog_burnt = fluid_model_local->GetTableLimitsProg().second; + for (unsigned long i_point = 0; i_point < nPointDomain; i_point++) { + for (unsigned long i_var = 0; i_var < nVar; i_var++) Solution[i_var] = 0.0; + + auto coords = geometry[i_mesh]->nodes->GetCoord(i_point); + + /* determine if our location 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 --- */ + } 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 --- */ + } else if ((point_loc > flame_thickness) && (point_loc <= flame_thickness + burnt_thickness)) { + scalar_init[I_PROGVAR] = prog_burnt; + n_points_burnt_local++; + + /* --- unburnt region downstream of the flame --- */ + } else { + scalar_init[I_PROGVAR] = prog_unburnt; + n_points_unburnt_local++; + } + + n_not_iterated_local += fluid_model_local->GetEnthFromTemp(&enth_inlet, prog_inlet, temp_inlet, enth_inlet); + scalar_init[I_ENTH] = enth_inlet; + + n_not_in_domain_local += fluid_model_local->GetLookUpTable()->LookUp_XY( + look_up_tags, look_up_data, scalar_init[I_PROGVAR], scalar_init[I_ENTH]); + + /* --- initialize the auxiliary transported scalars (not controlling variables) --- */ + for (int i_scalar = config->GetNControlVars(); i_scalar < config->GetNScalars(); ++i_scalar) { + scalar_init[i_scalar] = config->GetSpecies_Init()[i_scalar]; + } + + solver_container[i_mesh][SPECIES_SOL]->GetNodes()->SetSolution(i_point, scalar_init); + } + + solver_container[i_mesh][SPECIES_SOL]->InitiateComms(geometry[i_mesh], config, SOLUTION); + solver_container[i_mesh][SPECIES_SOL]->CompleteComms(geometry[i_mesh], config, SOLUTION); + + solver_container[i_mesh][FLOW_SOL]->InitiateComms(geometry[i_mesh], config, SOLUTION); + solver_container[i_mesh][FLOW_SOL]->CompleteComms(geometry[i_mesh], config, SOLUTION); + + solver_container[i_mesh][FLOW_SOL]->Preprocessing(geometry[i_mesh], solver_container[i_mesh], config, i_mesh, + NO_RK_ITER, RUNTIME_FLOW_SYS, false); + } + + /* --- sum up some counters over processes --- */ + SU2_MPI::Reduce(&n_not_in_domain_local, &n_not_in_domain_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm()); + SU2_MPI::Reduce(&n_not_iterated_local, &n_not_iterated_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm()); + SU2_MPI::Reduce(&n_points_unburnt_local, &n_points_unburnt_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm()); + SU2_MPI::Reduce(&n_points_burnt_local, &n_points_burnt_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm()); + SU2_MPI::Reduce(&n_points_flame_local, &n_points_flame_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm()); + + 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 (rank == MASTER_NODE && (n_not_in_domain_global > 0 || n_not_iterated_global > 0)) cout << endl; + + + if (rank == MASTER_NODE && n_not_in_domain_global > 0) + cout << " !!! Initial condition: Number of points outside of table domain: " << n_not_in_domain_global << " !!!" << endl; + + if (rank == MASTER_NODE && n_not_iterated_global > 0) + cout << " !!! Initial condition: Number of points in which enthalpy could not be iterated: " << n_not_iterated_global + << " !!!" << endl; + + if (rank == MASTER_NODE && (n_not_in_domain_global > 0 || n_not_iterated_global > 0)) cout << endl; + delete [] scalar_init; + } +} + +void CSpeciesFlameletSolver::SetPreconditioner(CGeometry* geometry, CSolver** solver_container, CConfig* config) { + unsigned short iVar; + unsigned long iPoint, total_index; + + su2double BetaInc2, Density, dRhodT, dRhodC, Temperature, Delta; + + bool variable_density = (config->GetKind_DensityModel() == INC_DENSITYMODEL::VARIABLE); + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + /*--- Access the primitive variables at this node. ---*/ + + Density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + BetaInc2 = solver_container[FLOW_SOL]->GetNodes()->GetBetaInc2(iPoint); + Temperature = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint); + + unsigned short nVar_Flow = solver_container[FLOW_SOL]->GetnVar(); + + su2double SolP = solver_container[FLOW_SOL]->LinSysSol[iPoint * nVar_Flow + 0]; + su2double SolT = solver_container[FLOW_SOL]->LinSysSol[iPoint * nVar_Flow + nDim + 1]; + + /*--- We need the derivative of the equation of state to build the + preconditioning matrix. For now, the only option is the ideal gas + law, but in the future, dRhodT should be in the fluid model. ---*/ + + if (variable_density) { + dRhodT = -Density / Temperature; + } else { + dRhodT = 0.0; + } + + /*--- Passive scalars have no impact on the density. ---*/ + + dRhodC = 0.0; + + /*--- Modify matrix diagonal with term including volume and time step. ---*/ + + su2double Vol = geometry->nodes->GetVolume(iPoint); + Delta = Vol / (config->GetCFLRedCoeff_Species() * solver_container[FLOW_SOL]->GetNodes()->GetDelta_Time(iPoint)); + + /*--- Calculating the inverse of the preconditioning matrix + that multiplies the time derivative during time integration. ---*/ + + if (implicit) { + for (iVar = 0; iVar < nVar; iVar++) { + total_index = iPoint * nVar + iVar; + + su2double scalar = nodes->GetSolution(iPoint, iVar); + + /*--- Compute the lag terms for the decoupled linear system from + the mean flow equations and add to the residual for the scalar. + In short, we are effectively making these terms explicit. ---*/ + + su2double artcompc1 = SolP * scalar / (Density * BetaInc2); + su2double artcompc2 = SolT * dRhodT * scalar / (Density); + + LinSysRes[total_index] += artcompc1 + artcompc2; + + /*--- Add the extra Jacobian term to the scalar system. ---*/ + + su2double Jaccomp = scalar * dRhodC + Density; + su2double JacTerm = Jaccomp * Delta; + + Jacobian.AddVal2Diag(iPoint, JacTerm); + } + } + } +} + +void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solver_container, + CNumerics** numerics_container, CConfig* config, unsigned short iMesh) { + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + bool axisymmetric = config->GetAxisymmetric(); + + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + + auto* first_numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num() * MAX_TERMS]; + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (auto i_point = 0u; i_point < nPointDomain; i_point++) { + /*--- Set primitive variables w/o reconstruction ---*/ + + first_numerics->SetPrimitive(flowNodes->GetPrimitive(i_point), nullptr); + + /*--- Set scalar variables w/o reconstruction ---*/ + + first_numerics->SetScalarVar(nodes->GetSolution(i_point), nullptr); + + first_numerics->SetDiffusionCoeff(nodes->GetDiffusivity(i_point), nodes->GetDiffusivity(i_point)); + + /*--- Set volume of the dual cell. ---*/ + + first_numerics->SetVolume(geometry->nodes->GetVolume(i_point)); + + /*--- Update scalar sources in the fluidmodel ---*/ + + /*--- Axisymmetry source term for the scalar equation. ---*/ + if (axisymmetric) { + /*--- Set y coordinate ---*/ + first_numerics->SetCoord(geometry->nodes->GetCoord(i_point), geometry->nodes->GetCoord(i_point)); + /*-- gradients necessary for axisymmetric flows only? ---*/ + first_numerics->SetScalarVarGradient(nodes->GetGradient(i_point), nullptr); + } + + /*--- Retrieve scalar sources from CVariable class and update numerics class data. ---*/ + first_numerics->SetScalarSources(nodes->GetScalarSources(i_point)); + + auto residual = first_numerics->ComputeResidual(config); + + /*--- Add Residual ---*/ + + LinSysRes.SubtractBlock(i_point, residual); + + /*--- Implicit part ---*/ + + if (implicit) Jacobian.SubtractBlock2Diag(i_point, residual.jacobian_i); + } + END_SU2_OMP_FOR +} + +void CSpeciesFlameletSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, + CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) { + string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + + su2double enth_inlet; + su2double temp_inlet = config->GetInlet_Ttotal(Marker_Tag); + const su2double* inlet_scalar_original = config->GetInlet_SpeciesVal(Marker_Tag); + su2double* inlet_scalar = const_cast(inlet_scalar_original); + + CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); + + /*--- We compute inlet enthalpy from the temperature and progress variable ---*/ + enth_inlet = inlet_scalar_original[I_ENTH]; + fluid_model_local->GetEnthFromTemp(&enth_inlet, inlet_scalar[I_PROGVAR], temp_inlet, inlet_scalar_original[I_ENTH]); + inlet_scalar[I_ENTH] = enth_inlet; + + /*--- Loop over all the vertices on this boundary marker ---*/ + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ + + if (!geometry->nodes->GetDomain(iPoint)) continue; + + if (config->GetSpecies_StrongBC()) { + nodes->SetSolution_Old(iPoint, inlet_scalar); + + LinSysRes.SetBlock_Zero(iPoint); + + for (auto iVar = 0; iVar < nVar; iVar++) { + nodes->SetVal_ResTruncError_Zero(iPoint, iVar); + } + + /*--- Includes 1 in the diagonal ---*/ + for (auto iVar = 0u; iVar < nVar; iVar++) { + auto total_index = iPoint * nVar + iVar; + Jacobian.DeleteValsRowi(total_index); + } + } else { + /*--- Normal vector for this vertex (negate for outward convention) ---*/ + + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); + + /*--- Allocate the value at the inlet ---*/ + + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + + /*--- Retrieve solution at the farfield boundary node ---*/ + + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + + /*--- Set various quantities in the solver class ---*/ + + conv_numerics->SetPrimitive(V_domain, V_inlet); + + /*--- Set the species variable state at the inlet. ---*/ + + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Inlet_SpeciesVars[val_marker][iVertex]); + + /*--- Set various other quantities in the solver class ---*/ + + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + + /*--- Compute the residual using an upwind scheme ---*/ + + auto residual = conv_numerics->ComputeResidual(config); + LinSysRes.AddBlock(iPoint, residual); + + /*--- Jacobian contribution for implicit integration ---*/ + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } + } + END_SU2_OMP_FOR +} + +void CSpeciesFlameletSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, + CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) { + /*--- Loop over all the vertices on this boundary marker ---*/ + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + /* strong zero flux Neumann boundary condition at the outlet */ + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ + + if (geometry->nodes->GetDomain(iPoint)) { + /*--- Allocate the value at the outlet ---*/ + + auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + for (auto iVar = 0u; iVar < nVar; iVar++) Solution[iVar] = nodes->GetSolution(Point_Normal, iVar); + + nodes->SetSolution_Old(iPoint, Solution); + + LinSysRes.SetBlock_Zero(iPoint); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + nodes->SetVal_ResTruncError_Zero(iPoint, iVar); + } + + /*--- Includes 1 in the diagonal ---*/ + for (auto iVar = 0u; iVar < nVar; iVar++) { + auto total_index = iPoint * nVar + iVar; + Jacobian.DeleteValsRowi(total_index); + } + } + } + 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) {} + +void CSpeciesFlameletSolver::BC_Isothermal_Wall(CGeometry* geometry, CSolver** solver_container, + CNumerics* conv_numerics, CNumerics* visc_numerics, CConfig* config, + unsigned short val_marker) { + unsigned long iVertex, iPoint, total_index; + + bool implicit = config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT; + string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + su2double temp_wall = config->GetIsothermal_Temperature(Marker_Tag); + CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); + CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes(); + + su2double enth_init, enth_wall, prog_wall; + unsigned long n_not_iterated = 0; + + /*--- Loop over all the vertices on this boundary marker ---*/ + + for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ + + if (geometry->nodes->GetDomain(iPoint)) { + /*--- Set enthalpy on the wall ---*/ + + prog_wall = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint)[I_PROGVAR]; + if(config->GetSpecies_StrongBC()){ + + /*--- Initial guess for enthalpy value ---*/ + + enth_init = nodes->GetSolution(iPoint, I_ENTH); + enth_wall = enth_init; + + n_not_iterated += fluid_model_local->GetEnthFromTemp(&enth_wall, prog_wall, temp_wall, enth_init); + + /*--- Impose the value of the enthalpy as a strong boundary + 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); + + LinSysRes(iPoint, I_ENTH) = 0.0; + + nodes->SetVal_ResTruncError_Zero(iPoint, I_ENTH); + + if (implicit) { + total_index = iPoint * nVar + I_ENTH; + + Jacobian.DeleteValsRowi(total_index); + } + } else { + /*--- Weak BC formulation ---*/ + const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); + + const su2double Area = GeometryToolbox::Norm(nDim, Normal); + + + const auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + + /*--- Get coordinates of i & nearest normal and compute distance ---*/ + + const auto Coord_i = geometry->nodes->GetCoord(iPoint); + const auto Coord_j = geometry->nodes->GetCoord(Point_Normal); + su2double Edge_Vector[MAXNDIM]; + GeometryToolbox::Distance(nDim, Coord_j, Coord_i, Edge_Vector); + su2double dist_ij_2 = GeometryToolbox::SquaredNorm(nDim, Edge_Vector); + su2double dist_ij = sqrt(dist_ij_2); + + /*--- Compute the normal gradient in temperature using Twall ---*/ + + su2double dTdn = -(flowNodes->GetTemperature(Point_Normal) - temp_wall)/dist_ij; + + /*--- Get thermal conductivity ---*/ + + su2double thermal_conductivity = flowNodes->GetThermalConductivity(iPoint); + + /*--- 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; + } + + } + } + if (rank == MASTER_NODE && n_not_iterated > 0) { + cout << " !!! Isothermal wall bc (" << Marker_Tag + << "): Number of points in which enthalpy could not be iterated: " << n_not_iterated << " !!!" << endl; + } +} + +void CSpeciesFlameletSolver::BC_ConjugateHeat_Interface(CGeometry* geometry, CSolver** solver_container, + CNumerics* conv_numerics, CConfig* config, + unsigned short val_marker) { + unsigned long iVertex, iPoint, total_index; + + bool implicit = config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT; + string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + su2double temp_wall = config->GetIsothermal_Temperature(Marker_Tag); + CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); + CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes(); + su2double enth_wall, enth_init, prog_wall; + unsigned long n_not_iterated = 0; + + /*--- Loop over all the vertices on this boundary marker ---*/ + + for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + temp_wall = GetConjugateHeatVariable(val_marker, iVertex, 0); + + /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ + + if (geometry->nodes->GetDomain(iPoint)) { + + if(config->GetSpecies_StrongBC()){ + + /*--- Initial guess for enthalpy ---*/ + + enth_init = nodes->GetSolution(iPoint, I_ENTH); + enth_wall = enth_init; + + /*--- Set enthalpy on the wall ---*/ + + prog_wall = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint)[I_PROGVAR]; + n_not_iterated += fluid_model_local->GetEnthFromTemp(&enth_wall, prog_wall, temp_wall, enth_init); + + /*--- Impose the value of the enthalpy as a strong boundary + 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); + + LinSysRes(iPoint, I_ENTH) = 0.0; + + nodes->SetVal_ResTruncError_Zero(iPoint, I_ENTH); + + if (implicit) { + total_index = iPoint * nVar + I_ENTH; + + Jacobian.DeleteValsRowi(total_index); + } + }else{ + /*--- Weak BC formulation ---*/ + const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); + + const su2double Area = GeometryToolbox::Norm(nDim, Normal); + + + const auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + + /*--- Get coordinates of i & nearest normal and compute distance ---*/ + + const auto Coord_i = geometry->nodes->GetCoord(iPoint); + const auto Coord_j = geometry->nodes->GetCoord(Point_Normal); + su2double Edge_Vector[MAXNDIM]; + GeometryToolbox::Distance(nDim, Coord_j, Coord_i, Edge_Vector); + su2double dist_ij_2 = GeometryToolbox::SquaredNorm(nDim, Edge_Vector); + su2double dist_ij = sqrt(dist_ij_2); + + /*--- Compute the normal gradient in temperature using Twall ---*/ + + su2double dTdn = -(flowNodes->GetTemperature(Point_Normal) - temp_wall)/dist_ij; + + /*--- Get thermal conductivity ---*/ + + su2double thermal_conductivity = flowNodes->GetThermalConductivity(iPoint); + + /*--- 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; + } + + } + } + if (rank == MASTER_NODE && n_not_iterated > 0) { + cout << " !!! CHT interface (" << Marker_Tag + << "): Number of points in which enthalpy could not be iterated: " << n_not_iterated << " !!!" << endl; + } +} + diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index 690de424eee1..17b6d2b867b2 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -199,9 +199,12 @@ void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); const bool energy = config->GetEnergy_Equation(); + const bool flamelet = (config->GetKind_FluidModel() == FLUID_FLAMELET); const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); - if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; + /*--- for the flamelet model, the temperature is saved to file, but the energy equation is off ---*/ + + if (incompressible && ((!energy) && (!weakly_coupled_heat) && (!flamelet))) skipVars--; /*--- Load data from the restart into correct containers. ---*/ @@ -553,7 +556,7 @@ void CSpeciesSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta const bool axisymmetric = config->GetAxisymmetric(); if (axisymmetric) { - CNumerics *numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS]; + CNumerics *numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num()*MAX_TERMS]; SU2_OMP_FOR_DYN(omp_chunk_size) for (auto iPoint = 0u; iPoint < nPointDomain; iPoint++) { diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index d68c550d1850..ffd7152346b8 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -129,8 +129,9 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); const bool energy = config->GetEnergy_Equation(); const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); + const bool flamelet = (config->GetKind_FluidModel() == FLUID_FLAMELET); - if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; + if (incompressible && ((!energy) && (!weakly_coupled_heat) && (!flamelet))) skipVars--; /*--- Load data from the restart into correct containers. ---*/ diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index c000575c2187..7b8b4faa66a4 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -59,7 +59,7 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do /*--- Set the value of the temperature directly ---*/ su2double Temperature = Solution(iPoint, nDim+1); - const auto check_temp = SetTemperature(iPoint,Temperature); + auto check_temp = SetTemperature(iPoint,Temperature); /*--- Use the fluid model to compute the new value of density. Note that the thermodynamic pressure is constant and decoupled @@ -69,6 +69,11 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do FluidModel->SetTDState_T(Temperature, scalar); + /*--- flamelet block ---*/ + Solution(iPoint,nDim+1) = FluidModel->GetTemperature(); + Temperature = Solution(iPoint,nDim+1); + check_temp = SetTemperature(iPoint, Temperature); + /*--- Set the value of the density ---*/ const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity()); diff --git a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp new file mode 100644 index 000000000000..d67e53d01ad2 --- /dev/null +++ b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp @@ -0,0 +1,57 @@ +/*! + * \file CSpeciesFlameletVariable.cpp + * \brief Definition of the variable fields for the flamelet class. + * \author D. Mayer, T. Economon, N. Beishuizen + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/variables/CSpeciesFlameletVariable.hpp" + +CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf, unsigned long npoint, + unsigned long ndim, unsigned long nvar, const CConfig* config) + : CSpeciesVariable(species_inf, npoint, ndim, nvar, config) { + + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) + for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution(iPoint, iVar) = species_inf[iVar]; + + Solution_Old = Solution; + + /*--- Allocate and initialize solution for the dual time strategy ---*/ + bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || + (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)); + + if (dual_time) { + Solution_time_n = Solution; + Solution_time_n1 = Solution; + } + + /*--- Allocate residual structures ---*/ + + Res_TruncError.resize(nPoint, nVar) = su2double(0.0); + + /* Allocate space for the source and scalars for visualization */ + + source_scalar.resize(nPoint, config->GetNScalars()) = su2double(0.0); + lookup_scalar.resize(nPoint, config->GetNLookups()) = su2double(0.0); + inside_table.resize(nPoint) = 0; +} diff --git a/TestCases/flamelet/README.md b/TestCases/flamelet/README.md new file mode 100644 index 000000000000..2094bf1ed244 --- /dev/null +++ b/TestCases/flamelet/README.md @@ -0,0 +1,5 @@ +# Flamelet regression cases + +## laminar_premixed_flame + +A laminar premixed flame stabilized on an isothermal burner with a fixed wall temperature using a low resolution look-up table for all thermo-chemical quantities diff --git a/TestCases/flamelet/laminar_premixed_flame/README.md b/TestCases/flamelet/laminar_premixed_flame/README.md new file mode 100644 index 000000000000..1f32fb9eff77 --- /dev/null +++ b/TestCases/flamelet/laminar_premixed_flame/README.md @@ -0,0 +1,5 @@ +# Laminar premixed flame + +Make sure to download the mesh and fgm from the test cases repo! + +TODO: Detailed explanation and description of expected results \ No newline at end of file diff --git a/TestCases/flamelet/laminar_premixed_flame/config.cfg b/TestCases/flamelet/laminar_premixed_flame/config.cfg new file mode 100644 index 000000000000..d2a711722c2a --- /dev/null +++ b/TestCases/flamelet/laminar_premixed_flame/config.cfg @@ -0,0 +1,131 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Laminar premixed flame stabilized on isothermal burner % +% Author: Nijso Beishuizen, Daniel Mayer % +% Institution: Bosch Thermotechnology % +% Date: 10/14/2022 % +% File Version 7.4.0 "Falcon" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER = INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +MATH_PROBLEM= DIRECT +RESTART_SOL = NO +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= VARIABLE +INC_ENERGY_EQUATION = YES +INC_DENSITY_INIT= 1.00 +INC_VELOCITY_INIT= (0.0, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 300.0 +INC_NONDIM= DIMENSIONAL +% -------------------- FLUID MODEL --------------------------------------- % +% +FLUID_MODEL= FLUID_FLAMELET +VISCOSITY_MODEL= FLAMELET +CONDUCTIVITY_MODEL= FLAMELET +DIFFUSIVITY_MODEL= FLAMELET +KIND_SCALAR_MODEL= FLAMELET +FILENAME_LUT= fgm.drg +FLAME_OFFSET= 0 0.0035 0 +FLAME_THICKNESS= 0.2e-3 +FLAME_NORMAL= 0 1 0 +FLAME_BURNT_THICKNESS = 1 +PRANDTL_LAM= 1.0 +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +SPECIES_USE_STRONG_BC= NO +SPECIES_INIT = (0,-210000, 0, 0) + +% Passive reactants in flamelet problem +USER_SCALAR_NAMES = (Y-CO, Y-NOx) +USER_SOURCE_NAMES = ( \ + ProdRatePos_Y-CO, ProdRateNegScaled_Y-CO ,\ + ProdRateTot_Y-NOx, NULL +) + +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR +MUSCL_SPECIES= NO +SLOPE_LIMITER_SPECIES= NONE +TIME_DISCRE_SPECIES= EULER_IMPLICIT +% SCALAR CLIPPING +SPECIES_CLIPPING= NO +%SPECIES_CLIPPING_MIN= 0 -2400000 0 0 +%SPECIES_CLIPPING_MAX= 0.2500 400000 1 1 +MARKER_INLET_SPECIES = (inlet, 0, 0, 0, 0) +CFL_REDUCTION_SPECIES= 1.0 +LOOKUP_NAMES = (HeatRelease) +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= (wall, 500) +MARKER_SYM= (symmetry, symmetry_2-part-fluid) +INC_INLET_TYPE= VELOCITY_INLET +MARKER_INLET = (inlet, 300.000, 0.3, 0.0, 1.0, 0.0) +INC_OUTLET_TYPE= PRESSURE_OUTLET +INC_INLET_DAMPING = 0.05 +INC_OUTLET_DAMPING = 0.05 +MARKER_OUTLET= (outlet, 0.0) +MARKER_PLOTTING= ( wall ) +MARKER_MONITORING= ( wall ) +MARKER_ANALYZE= ( inlet,outlet ) +MARKER_ANALYZE_AVERAGE = AREA +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER =50.0 +CFL_ADAPT= NO +ITER=30000 +OUTPUT_WRT_FREQ= 500 +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-10 +LINEAR_SOLVER_ITER=20 +% +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 0 +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= NO +SLOPE_LIMITER_FLOW = NONE +VENKAT_LIMITER_COEFF= 10.0 +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -32 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 +SCREEN_OUTPUT = INNER_ITER WALL_TIME LINSOL_ITER LINSOL_RESIDUAL LINSOL_ITER_SPECIES LINSOL_RESIDUAL_SPECIES RMS_SPECIES RMS_PRESSURE RMS_TEMPERATURE RMS_VELOCITY_X RMS_PROGRESS_VARIABLE RMS_TOTAL_ENTHALPY N_TABLE_MISSES +HISTORY_OUTPUT = RMS_RES AERO_COEFF FLOW_COEFF FLOW_COEFF_SURF +VOLUME_OUTPUT = SOLUTION PRIMITIVE SOURCE RESIDUAL SENSITIVITY LOOKUP TIMESTEP +%FLAMELET_TABLE_OUTPUT = Conductivity ViscosityDyn Diffusivity Density +CONV_FIELD = RMS_PRESSURE +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= CGNS +MESH_FILENAME = mesh.cgns +SOLUTION_FILENAME= solution.dat +OUTPUT_FILES = (RESTART,PARAVIEW) +TABULAR_FORMAT = CSV +CONV_FILENAME= history +RESTART_FILENAME= restart.dat +VOLUME_FILENAME= flow +WRT_PERFORMANCE = YES +SCREEN_WRT_FREQ_INNER = 1 +SCREEN_WRT_FREQ_OUTER = 1 \ No newline at end of file diff --git a/config_template.cfg b/config_template.cfg index 0463ce3357f7..058225f9406c 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -317,7 +317,7 @@ SEMI_SPAN= 0.0 % ---- NONEQUILIBRIUM GAS, IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS, CoolProp library -------% % % Fluid model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS, -% CONSTANT_DENSITY, INC_IDEAL_GAS, INC_IDEAL_GAS_POLY, MUTATIONPP, SU2_NONEQ, FLUID_MIXTURE, COOLPROP) +% CONSTANT_DENSITY, INC_IDEAL_GAS, INC_IDEAL_GAS_POLY, MUTATIONPP, SU2_NONEQ, FLUID_MIXTURE, COOLPROP, FLUID_FLAMELET) FLUID_MODEL= STANDARD_AIR % To find all available fluid name for CoolProp library, clikc the following link: % http://www.coolprop.org/fluid_properties/PurePseudoPure.html#list-of-fluids @@ -706,7 +706,7 @@ TIME_DISCRE_RADIATION = EULER_IMPLICIT % --------------------- SPECIES TRANSPORT SIMULATION --------------------------% % -% Specify scalar transport model (NONE, SPECIES_TRANSPORT) +% Specify scalar transport model (NONE, SPECIES_TRANSPORT, FLAMELET) KIND_SCALAR_MODEL= NONE % % Mass diffusivity model (CONSTANT_DIFFUSIVITY)