diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 4fcfea1531d6..1e501b3e4ea9 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -531,28 +531,34 @@ class CConfig { Kind_ConvNumScheme_AdjTurb, /*!< \brief Centered or upwind scheme for the adjoint turbulence model. */ Kind_ConvNumScheme_Species, /*!< \brief Centered or upwind scheme for the species model. */ Kind_ConvNumScheme_Template, /*!< \brief Centered or upwind scheme for the level set equation. */ + Kind_FEM, /*!< \brief Finite element scheme for the flow equations. */ + Kind_FEM_Flow, /*!< \brief Finite element scheme for the flow equations. */ + Kind_Matrix_Coloring; /*!< \brief Type of matrix coloring for sparse Jacobian computation. */ + + CENTERED Kind_Centered, /*!< \brief Centered scheme. */ Kind_Centered_Flow, /*!< \brief Centered scheme for the flow equations. */ Kind_Centered_AdjFlow, /*!< \brief Centered scheme for the adjoint flow equations. */ Kind_Centered_Turb, /*!< \brief Centered scheme for the turbulence model. */ Kind_Centered_AdjTurb, /*!< \brief Centered scheme for the adjoint turbulence model. */ Kind_Centered_Species, /*!< \brief Centered scheme for the species model. */ - Kind_Centered_Template, /*!< \brief Centered scheme for the template model. */ + Kind_Centered_Template; /*!< \brief Centered scheme for the template model. */ + + + FEM_SHOCK_CAPTURING_DG Kind_FEM_Shock_Capturing_DG; /*!< \brief Shock capturing method for the FEM DG solver. */ + BGS_RELAXATION Kind_BGS_RelaxMethod; /*!< \brief Kind of relaxation method for Block Gauss Seidel method in FSI problems. */ + bool ReconstructionGradientRequired; /*!< \brief Enable or disable a second gradient calculation for upwind reconstruction only. */ + bool LeastSquaresRequired; /*!< \brief Enable or disable memory allocation for least-squares gradient methods. */ + bool Energy_Equation; /*!< \brief Solve the energy equation for incompressible flows. */ + + UPWIND Kind_Upwind, /*!< \brief Upwind scheme. */ Kind_Upwind_Flow, /*!< \brief Upwind scheme for the flow equations. */ Kind_Upwind_AdjFlow, /*!< \brief Upwind scheme for the adjoint flow equations. */ Kind_Upwind_Turb, /*!< \brief Upwind scheme for the turbulence model. */ Kind_Upwind_AdjTurb, /*!< \brief Upwind scheme for the adjoint turbulence model. */ Kind_Upwind_Species, /*!< \brief Upwind scheme for the species model. */ - Kind_Upwind_Template, /*!< \brief Upwind scheme for the template model. */ - Kind_FEM, /*!< \brief Finite element scheme for the flow equations. */ - Kind_FEM_Flow, /*!< \brief Finite element scheme for the flow equations. */ - Kind_Matrix_Coloring; /*!< \brief Type of matrix coloring for sparse Jacobian computation. */ - FEM_SHOCK_CAPTURING_DG Kind_FEM_Shock_Capturing_DG; /*!< \brief Shock capturing method for the FEM DG solver. */ - BGS_RELAXATION Kind_BGS_RelaxMethod; /*!< \brief Kind of relaxation method for Block Gauss Seidel method in FSI problems. */ - bool ReconstructionGradientRequired; /*!< \brief Enable or disable a second gradient calculation for upwind reconstruction only. */ - bool LeastSquaresRequired; /*!< \brief Enable or disable memory allocation for least-squares gradient methods. */ - bool Energy_Equation; /*!< \brief Solve the energy equation for incompressible flows. */ + Kind_Upwind_Template; /*!< \brief Upwind scheme for the template model. */ bool MUSCL, /*!< \brief MUSCL scheme .*/ MUSCL_Flow, /*!< \brief MUSCL scheme for the flow equations.*/ @@ -1103,7 +1109,7 @@ class CConfig { hs_axes[3], /*!< \brief principal axes (x, y, z) of the ellipsoid containing the heat source. */ hs_center[3]; /*!< \brief position of the center of the heat source. */ - unsigned short Riemann_Solver_FEM; /*!< \brief Riemann solver chosen for the DG method. */ + UPWIND Riemann_Solver_FEM; /*!< \brief Riemann solver chosen for the DG method. */ su2double Quadrature_Factor_Straight; /*!< \brief Factor applied during quadrature of elements with a constant Jacobian. */ su2double Quadrature_Factor_Curved; /*!< \brief Factor applied during quadrature of elements with a non-constant Jacobian. */ su2double Quadrature_Factor_Time_ADER_DG; /*!< \brief Factor applied during quadrature in time for ADER-DG. */ @@ -1274,7 +1280,7 @@ class CConfig { void addStringListOption(const string name, unsigned short & num_marker, string* & option_field); - void addConvectOption(const string name, unsigned short & space_field, unsigned short & centered_field, unsigned short & upwind_field); + void addConvectOption(const string name, unsigned short & space_field, CENTERED & centered_field, UPWIND & upwind_field); void addConvectFEMOption(const string name, unsigned short & space_field, unsigned short & fem_field); @@ -2300,8 +2306,8 @@ class CConfig { * \param[in] val_muscl - Define if we apply a MUSCL scheme or not. * \param[in] val_kind_fem - If FEM, what kind of FEM discretization. */ - void SetKind_ConvNumScheme(unsigned short val_kind_convnumscheme, unsigned short val_kind_centered, - unsigned short val_kind_upwind, LIMITER val_kind_slopelimit, + void SetKind_ConvNumScheme(unsigned short val_kind_convnumscheme, CENTERED val_kind_centered, + UPWIND val_kind_upwind, LIMITER val_kind_slopelimit, bool val_muscl, unsigned short val_kind_fem); /*! @@ -3688,7 +3694,7 @@ class CConfig { */ bool GetAUSMMethod(void) const { switch (Kind_Upwind_Flow) { - case AUSM : case AUSMPLUSUP: case AUSMPLUSUP2: case AUSMPWPLUS: + case UPWIND::AUSM : case UPWIND::AUSMPLUSUP: case UPWIND::AUSMPLUSUP2: case UPWIND::AUSMPWPLUS: return true; default: return false; @@ -4341,7 +4347,7 @@ class CConfig { * linearized) that is being solved. * \return Kind of center scheme for the convective terms. */ - unsigned short GetKind_Centered(void) const { return Kind_Centered; } + CENTERED GetKind_Centered(void) const { return Kind_Centered; } /*! * \brief Get kind of upwind scheme for the convective terms. @@ -4350,7 +4356,7 @@ class CConfig { * linearized) that is being solved. * \return Kind of upwind scheme for the convective terms. */ - unsigned short GetKind_Upwind(void) const { return Kind_Upwind; } + UPWIND GetKind_Upwind(void) const { return Kind_Upwind; } /*! * \brief Get if the upwind scheme used MUSCL or not. @@ -4521,7 +4527,7 @@ class CConfig { * during the computation. * \return Kind of center convective numerical scheme for the flow equations. */ - ENUM_CENTERED GetKind_Centered_Flow(void) const { return static_cast(Kind_Centered_Flow); } + CENTERED GetKind_Centered_Flow(void) const { return Kind_Centered_Flow; } /*! * \brief Get the kind of center convective numerical scheme for the plasma equations. @@ -4537,7 +4543,7 @@ class CConfig { * during the computation. * \return Kind of upwind convective numerical scheme for the flow equations. */ - unsigned short GetKind_Upwind_Flow(void) const { return Kind_Upwind_Flow; } + UPWIND GetKind_Upwind_Flow(void) const { return Kind_Upwind_Flow; } /*! * \brief Get the kind of finite element convective numerical scheme for the flow equations. @@ -4666,7 +4672,7 @@ class CConfig { * during the computation. * \return Kind of center convective numerical scheme for the adjoint flow equations. */ - unsigned short GetKind_Centered_AdjFlow(void) const { return Kind_Centered_AdjFlow; } + CENTERED GetKind_Centered_AdjFlow(void) const { return Kind_Centered_AdjFlow; } /*! * \brief Get the kind of upwind convective numerical scheme for the adjoint flow equations. @@ -4674,7 +4680,7 @@ class CConfig { * during the computation. * \return Kind of upwind convective numerical scheme for the adjoint flow equations. */ - unsigned short GetKind_Upwind_AdjFlow(void) const { return Kind_Upwind_AdjFlow; } + UPWIND GetKind_Upwind_AdjFlow(void) const { return Kind_Upwind_AdjFlow; } /*! * \brief Value of the calibrated constant for the high order method (center scheme). @@ -4718,7 +4724,7 @@ class CConfig { * during the computation. * \return Kind of center convective numerical scheme for the turbulence equations. */ - unsigned short GetKind_Centered_Turb(void) const { return Kind_Centered_Turb; } + CENTERED GetKind_Centered_Turb(void) const { return Kind_Centered_Turb; } /*! * \brief Get the kind of upwind convective numerical scheme for the turbulence equations. @@ -4726,7 +4732,7 @@ class CConfig { * during the computation. * \return Kind of upwind convective numerical scheme for the turbulence equations. */ - unsigned short GetKind_Upwind_Turb(void) const { return Kind_Upwind_Turb; } + UPWIND GetKind_Upwind_Turb(void) const { return Kind_Upwind_Turb; } /*! * \brief Get the kind of integration scheme (explicit or implicit) @@ -4770,7 +4776,7 @@ class CConfig { * during the computation. * \return Kind of center convective numerical scheme for the Species equations. */ - unsigned short GetKind_Centered_Species() const { return Kind_Centered_Species; } + CENTERED GetKind_Centered_Species() const { return Kind_Centered_Species; } /*! * \brief Get the kind of upwind convective numerical scheme for the Species equations. @@ -4778,7 +4784,22 @@ class CConfig { * during the computation. * \return Kind of upwind convective numerical scheme for the Species equations. */ - unsigned short GetKind_Upwind_Species() const { return Kind_Upwind_Species; } + UPWIND GetKind_Upwind_Species() const { return Kind_Upwind_Species; } + + /*! + * \brief Returns true if bounded scalar mode is on for species transport. + */ + bool GetBounded_Species() const { return (Kind_Upwind_Species == UPWIND::BOUNDED_SCALAR); } + + /*! + * \brief Returns true if bounded scalar mode is on for turbulence transport. + */ + bool GetBounded_Turb() const { return (Kind_Upwind_Turb == UPWIND::BOUNDED_SCALAR); } + + /*! + * \brief Returns true if bounded scalar mode is used for any equation. + */ + bool GetBounded_Scalar() const { return GetBounded_Species() || GetBounded_Turb(); } /*! * \brief Get the kind of convective numerical scheme for the heat equation. @@ -4794,7 +4815,7 @@ class CConfig { * during the computation. * \return Kind of center convective numerical scheme for the adjoint turbulence equations. */ - unsigned short GetKind_Centered_AdjTurb(void) const { return Kind_Centered_AdjTurb; } + CENTERED GetKind_Centered_AdjTurb(void) const { return Kind_Centered_AdjTurb; } /*! * \brief Get the kind of upwind convective numerical scheme for the adjoint turbulence equations. @@ -4802,7 +4823,7 @@ class CConfig { * during the computation. * \return Kind of upwind convective numerical scheme for the adjoint turbulence equations. */ - unsigned short GetKind_Upwind_AdjTurb(void) const { return Kind_Upwind_AdjTurb; } + UPWIND GetKind_Upwind_AdjTurb(void) const { return Kind_Upwind_AdjTurb; } /*! * \brief Provides information about the way in which the turbulence will be treated by the @@ -8898,7 +8919,7 @@ class CConfig { * \note This value is obtained from the config file, and it is constant during the computation. * \return Kind of Riemann solver for the DG method (FEM flow solver). */ - unsigned short GetRiemann_Solver_FEM(void) const { return Riemann_Solver_FEM; } + UPWIND GetRiemann_Solver_FEM(void) const { return Riemann_Solver_FEM; } /*! * \brief Get the factor applied during quadrature of straight elements. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 7d076416cd83..152a3829d27f 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -809,67 +809,69 @@ static const MapType Gust_Dir_Map = { /*! * \brief Types of centered spatial discretizations */ -enum ENUM_CENTERED { - NO_CENTERED = 0, /*!< \brief No centered scheme is used. */ - JST = 1, /*!< \brief Jameson-Smith-Turkel centered numerical method. */ - LAX = 2, /*!< \brief Lax-Friedrich centered numerical method. */ - JST_MAT = 3, /*!< \brief JST with matrix dissipation. */ - JST_KE = 4 /*!< \brief Kinetic Energy preserving Jameson-Smith-Turkel centered numerical method. */ +enum class CENTERED { + NONE, /*!< \brief No centered scheme is used. */ + JST, /*!< \brief Jameson-Smith-Turkel centered numerical method. */ + LAX, /*!< \brief Lax-Friedrich centered numerical method. */ + JST_MAT, /*!< \brief JST with matrix dissipation. */ + JST_KE /*!< \brief Kinetic Energy preserving Jameson-Smith-Turkel centered numerical method. */ }; -static const MapType Centered_Map = { - MakePair("NONE", NO_CENTERED) - MakePair("JST", JST) - MakePair("JST_KE", JST_KE) - MakePair("JST_MAT", JST_MAT) - MakePair("LAX-FRIEDRICH", LAX) +static const MapType Centered_Map = { + MakePair("NONE", CENTERED::NONE) + MakePair("JST", CENTERED::JST) + MakePair("JST_KE", CENTERED::JST_KE) + MakePair("JST_MAT", CENTERED::JST_MAT) + MakePair("LAX-FRIEDRICH", CENTERED::LAX) }; -// If you add to ENUM_UPWIND, you must also add the option to ENUM_CONVECTIVE +// If you add to UPWIND, you must also add the option to ENUM_CONVECTIVE /*! * \brief Types of upwind spatial discretizations */ -enum ENUM_UPWIND { - NO_UPWIND = 0, /*!< \brief No upwind scheme is used. */ - ROE = 1, /*!< \brief Roe's upwind numerical method. */ - SCALAR_UPWIND = 2, /*!< \brief Scalar upwind numerical method. */ - AUSM = 3, /*!< \brief AUSM numerical method. */ - HLLC = 4, /*!< \brief HLLC numerical method. */ - SW = 5, /*!< \brief Steger-Warming method. */ - MSW = 6, /*!< \brief Modified Steger-Warming method. */ - TURKEL = 7, /*!< \brief Roe-Turkel's upwind numerical method. */ - SLAU = 8, /*!< \brief Simple Low-Dissipation AUSM numerical method. */ - CUSP = 9, /*!< \brief Convective upwind and split pressure numerical method. */ - CONVECTIVE_TEMPLATE = 10, /*!< \brief Template for new numerical method . */ - L2ROE = 11, /*!< \brief L2ROE numerical method . */ - LMROE = 12, /*!< \brief Rieper's Low Mach ROE numerical method . */ - SLAU2 = 13, /*!< \brief Simple Low-Dissipation AUSM 2 numerical method. */ - FDS = 14, /*!< \brief Flux difference splitting upwind method (incompressible flows). */ - LAX_FRIEDRICH = 15, /*!< \brief Lax-Friedrich numerical method. */ - AUSMPLUSUP = 16, /*!< \brief AUSM+ -up numerical method (All Speed) */ - AUSMPLUSUP2 = 17, /*!< \brief AUSM+ -up2 numerical method (All Speed) */ - AUSMPWPLUS = 18 /*!< \brief AUSMplus numerical method. (MAYBE for TNE2 ONLY)*/ -}; -static const MapType Upwind_Map = { - MakePair("NONE", NO_UPWIND) - MakePair("ROE", ROE) - MakePair("TURKEL_PREC", TURKEL) - MakePair("AUSM", AUSM) - MakePair("AUSMPLUSUP", AUSMPLUSUP) - MakePair("AUSMPLUSUP2", AUSMPLUSUP2) - MakePair("AUSMPWPLUS", AUSMPWPLUS) - MakePair("SLAU", SLAU) - MakePair("HLLC", HLLC) - MakePair("SW", SW) - MakePair("MSW", MSW) - MakePair("CUSP", CUSP) - MakePair("SCALAR_UPWIND", SCALAR_UPWIND) - MakePair("CONVECTIVE_TEMPLATE", CONVECTIVE_TEMPLATE) - MakePair("L2ROE", L2ROE) - MakePair("LMROE", LMROE) - MakePair("SLAU2", SLAU2) - MakePair("FDS", FDS) - MakePair("LAX-FRIEDRICH", LAX_FRIEDRICH) +enum class UPWIND { + NONE, /*!< \brief No upwind scheme is used. */ + ROE, /*!< \brief Roe's upwind numerical method. */ + SCALAR_UPWIND, /*!< \brief Scalar upwind numerical method. */ + AUSM, /*!< \brief AUSM numerical method. */ + HLLC, /*!< \brief HLLC numerical method. */ + SW, /*!< \brief Steger-Warming method. */ + MSW, /*!< \brief Modified Steger-Warming method. */ + TURKEL, /*!< \brief Roe-Turkel's upwind numerical method. */ + SLAU, /*!< \brief Simple Low-Dissipation AUSM numerical method. */ + CUSP, /*!< \brief Convective upwind and split pressure numerical method. */ + CONVECTIVE_TEMPLATE, /*!< \brief Template for new numerical method . */ + L2ROE, /*!< \brief L2ROE numerical method . */ + LMROE, /*!< \brief Rieper's Low Mach ROE numerical method . */ + SLAU2, /*!< \brief Simple Low-Dissipation AUSM 2 numerical method. */ + FDS, /*!< \brief Flux difference splitting upwind method (incompressible flows). */ + LAX_FRIEDRICH, /*!< \brief Lax-Friedrich numerical method. */ + AUSMPLUSUP, /*!< \brief AUSM+ -up numerical method (All Speed) */ + AUSMPLUSUP2, /*!< \brief AUSM+ -up2 numerical method (All Speed) */ + AUSMPWPLUS, /*!< \brief AUSMplus numerical method. (MAYBE for TNE2 ONLY)*/ + BOUNDED_SCALAR /*!< \brief Scalar advection numerical method. */ +}; +static const MapType Upwind_Map = { + MakePair("NONE", UPWIND::NONE) + MakePair("ROE", UPWIND::ROE) + MakePair("TURKEL_PREC", UPWIND::TURKEL) + MakePair("AUSM", UPWIND::AUSM) + MakePair("AUSMPLUSUP", UPWIND::AUSMPLUSUP) + MakePair("AUSMPLUSUP2", UPWIND::AUSMPLUSUP2) + MakePair("AUSMPWPLUS", UPWIND::AUSMPWPLUS) + MakePair("SLAU", UPWIND::SLAU) + MakePair("HLLC", UPWIND::HLLC) + MakePair("SW", UPWIND::SW) + MakePair("MSW", UPWIND::MSW) + MakePair("CUSP", UPWIND::CUSP) + MakePair("SCALAR_UPWIND", UPWIND::SCALAR_UPWIND) + MakePair("BOUNDED_SCALAR", UPWIND::BOUNDED_SCALAR) + MakePair("CONVECTIVE_TEMPLATE", UPWIND::CONVECTIVE_TEMPLATE) + MakePair("L2ROE", UPWIND::L2ROE) + MakePair("LMROE", UPWIND::LMROE) + MakePair("SLAU2", UPWIND::SLAU2) + MakePair("FDS", UPWIND::FDS) + MakePair("LAX-FRIEDRICH", UPWIND::LAX_FRIEDRICH) }; /*! diff --git a/Common/include/option_structure.inl b/Common/include/option_structure.inl index 4973de92cb36..5e09efe9aa24 100644 --- a/Common/include/option_structure.inl +++ b/Common/include/option_structure.inl @@ -406,11 +406,11 @@ public: class COptionConvect : public COptionBase { string name; // identifier for the option unsigned short & space; - unsigned short & centered; - unsigned short & upwind; + CENTERED & centered; + UPWIND & upwind; public: - COptionConvect(string option_field_name, unsigned short & space_field, unsigned short & centered_field, unsigned short & upwind_field) + COptionConvect(string option_field_name, unsigned short & space_field, CENTERED & centered_field, UPWIND & upwind_field) : name(option_field_name), space(space_field), centered(centered_field), upwind(upwind_field) { } string SetValue(const vector& option_value) override { @@ -424,13 +424,13 @@ public: if (Centered_Map.count(option_value[0])) { this->space = Space_Map.find("SPACE_CENTERED")->second; this->centered = Centered_Map.find(option_value[0])->second; - this->upwind = NO_UPWIND; + this->upwind = UPWIND::NONE; return ""; } if (Upwind_Map.count(option_value[0])) { this->space = Space_Map.find("SPACE_UPWIND")->second; this->upwind = Upwind_Map.find(option_value[0])->second; - this->centered = NO_CENTERED; + this->centered = CENTERED::NONE; return ""; } // Make them defined in case something weird happens @@ -440,8 +440,8 @@ public: } void SetDefault() override { - this->centered = NO_CENTERED; - this->upwind = NO_UPWIND; + this->centered = CENTERED::NONE; + this->upwind = UPWIND::NONE; this->space = NO_CONVECTIVE; } }; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 61795b75226d..c55860989809 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -414,7 +414,7 @@ void CConfig::addStringListOption(const string name, unsigned short & num_marker option_map.insert(pair(name, val)); } -void CConfig::addConvectOption(const string name, unsigned short & space_field, unsigned short & centered_field, unsigned short & upwind_field) { +void CConfig::addConvectOption(const string name, unsigned short & space_field, CENTERED & centered_field, UPWIND & upwind_field) { assert(option_map.find(name) == option_map.end()); all_options.insert(pair(name, true)); COptionBase* val = new COptionConvect(name, space_field, centered_field, upwind_field); @@ -2311,7 +2311,7 @@ void CConfig::SetConfig_Options() { /*--- Options related to the finite element flow solver---*/ /* DESCRIPTION: Riemann solver used for DG (ROE, LAX-FRIEDRICH, AUSM, AUSMPW+, HLLC, VAN_LEER) */ - addEnumOption("RIEMANN_SOLVER_FEM", Riemann_Solver_FEM, Upwind_Map, ROE); + addEnumOption("RIEMANN_SOLVER_FEM", Riemann_Solver_FEM, Upwind_Map, UPWIND::ROE); /* DESCRIPTION: Constant factor applied for quadrature with straight elements (2.0 by default) */ addDoubleOption("QUADRATURE_FACTOR_STRAIGHT_FEM", Quadrature_Factor_Straight, 2.0); /* DESCRIPTION: Constant factor applied for quadrature with curved elements (3.0 by default) */ @@ -3897,18 +3897,16 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i CURRENT_FUNCTION); } - if (!ideal_gas && !nemo) { - if (Kind_Upwind_Flow != ROE && Kind_Upwind_Flow != HLLC && Kind_Centered_Flow != JST) { - SU2_MPI::Error( - "Only ROE Upwind, HLLC Upwind scheme, and JST scheme can be used for Non-Ideal Compressible Fluids", - CURRENT_FUNCTION); - } + if (!ideal_gas && !nemo) { + if (Kind_Upwind_Flow != UPWIND::ROE && Kind_Upwind_Flow != UPWIND::HLLC && Kind_Centered_Flow != CENTERED::JST) { + SU2_MPI::Error("Only ROE Upwind, HLLC Upwind scheme, and JST scheme can be used for Non-Ideal Compressible Fluids", CURRENT_FUNCTION); } + } - if (nemo) { - if (Kind_Upwind_Flow == AUSMPWPLUS) - SU2_MPI::Error("AUSMPW+ is extremely unstable. Feel free to fix me!", CURRENT_FUNCTION); - } + if (nemo){ + if (Kind_Upwind_Flow == UPWIND::AUSMPWPLUS) + SU2_MPI::Error("AUSMPW+ is extremely unstable. Feel free to fix me!", CURRENT_FUNCTION); + } if (GetBoolTurbomachinery()) { nBlades = new su2double[nZone]; @@ -5433,6 +5431,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } } // species transport checks + + if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) { + SU2_MPI::Error("BOUNDED_SCALAR discretization can only be used for incompressible problems.", CURRENT_FUNCTION); + } } void CConfig::SetMarkers(SU2_COMPONENT val_software) { @@ -6594,7 +6596,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { (Kind_Solver == MAIN_SOLVER::DISC_ADJ_EULER) || (Kind_Solver == MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES) || (Kind_Solver == MAIN_SOLVER::DISC_ADJ_RANS) ) { if (Kind_ConvNumScheme_Flow == SPACE_CENTERED) { - if (Kind_Centered_Flow == LAX) { + if (Kind_Centered_Flow == CENTERED::LAX) { cout << "Lax-Friedrich scheme (1st order in space) for the flow inviscid terms.\n"; cout << "Lax viscous coefficients (1st): " << Kappa_1st_Flow << ".\n"; cout << "First order integration." << endl; @@ -6607,21 +6609,21 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { } if (Kind_ConvNumScheme_Flow == SPACE_UPWIND) { - if (Kind_Upwind_Flow == ROE) cout << "Roe (with entropy fix = "<< EntropyFix_Coeff <<") solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == TURKEL) cout << "Roe-Turkel solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == AUSM) cout << "AUSM solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == HLLC) cout << "HLLC solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == SW) cout << "Steger-Warming solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == MSW) cout << "Modified Steger-Warming solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == CUSP) cout << "CUSP solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == L2ROE) cout << "L2ROE Low Mach ROE solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == LMROE) cout << "Rieper Low Mach ROE solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == SLAU) cout << "Simple Low-Dissipation AUSM solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == SLAU2) cout << "Simple Low-Dissipation AUSM 2 solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == FDS) cout << "Flux difference splitting (FDS) upwind scheme for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == AUSMPLUSUP) cout << "AUSM+-up solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == AUSMPLUSUP2) cout << "AUSM+-up2 solver for the flow inviscid terms."<< endl; - if (Kind_Upwind_Flow == AUSMPWPLUS) cout << "AUSMPWPLUS solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::ROE) cout << "Roe (with entropy fix = "<< EntropyFix_Coeff <<") solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::TURKEL) cout << "Roe-Turkel solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::AUSM) cout << "AUSM solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::HLLC) cout << "HLLC solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::SW) cout << "Steger-Warming solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::MSW) cout << "Modified Steger-Warming solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::CUSP) cout << "CUSP solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::L2ROE) cout << "L2ROE Low Mach ROE solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::LMROE) cout << "Rieper Low Mach ROE solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::SLAU) cout << "Simple Low-Dissipation AUSM solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::SLAU2) cout << "Simple Low-Dissipation AUSM 2 solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::FDS) cout << "Flux difference splitting (FDS) upwind scheme for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::AUSMPLUSUP) cout << "AUSM+-up solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::AUSMPLUSUP2) cout << "AUSM+-up2 solver for the flow inviscid terms."<< endl; + if (Kind_Upwind_Flow == UPWIND::AUSMPWPLUS) cout << "AUSMPWPLUS solver for the flow inviscid terms."<< endl; if (Kind_Solver == MAIN_SOLVER::EULER || Kind_Solver == MAIN_SOLVER::DISC_ADJ_EULER || Kind_Solver == MAIN_SOLVER::NAVIER_STOKES || Kind_Solver == MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES || @@ -6647,7 +6649,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { if ((Kind_Solver == MAIN_SOLVER::RANS) || (Kind_Solver == MAIN_SOLVER::DISC_ADJ_RANS)) { if (Kind_ConvNumScheme_Turb == SPACE_UPWIND) { - if (Kind_Upwind_Turb == SCALAR_UPWIND) cout << "Scalar upwind solver for the turbulence model." << endl; + if (Kind_Upwind_Turb == UPWIND::SCALAR_UPWIND) cout << "Scalar upwind solver for the turbulence model." << endl; if (MUSCL_Turb) { PrintLimiterInfo(Kind_SlopeLimit_Turb); } else { @@ -6659,21 +6661,21 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { if ((Kind_Solver == MAIN_SOLVER::ADJ_EULER) || (Kind_Solver == MAIN_SOLVER::ADJ_NAVIER_STOKES) || (Kind_Solver == MAIN_SOLVER::ADJ_RANS)) { if (Kind_ConvNumScheme_AdjFlow == SPACE_CENTERED) { - if (Kind_Centered_AdjFlow == JST) { + if (Kind_Centered_AdjFlow == CENTERED::JST) { cout << "Jameson-Schmidt-Turkel scheme for the adjoint inviscid terms."<< endl; cout << "JST viscous coefficients (1st, 2nd, & 4th): " << Kappa_1st_AdjFlow << ", " << Kappa_2nd_AdjFlow << ", " << Kappa_4th_AdjFlow <<"."<< endl; cout << "The method includes a grid stretching correction (p = 0.3)."<< endl; cout << "Second order integration." << endl; } - if (Kind_Centered_AdjFlow == LAX) { + if (Kind_Centered_AdjFlow == CENTERED::LAX) { cout << "Lax-Friedrich scheme for the adjoint inviscid terms."<< endl; cout << "First order integration." << endl; } } if (Kind_ConvNumScheme_AdjFlow == SPACE_UPWIND) { - if (Kind_Upwind_AdjFlow == ROE) cout << "Roe (with entropy fix = "<< EntropyFix_Coeff <<") solver for the adjoint inviscid terms."<< endl; + if (Kind_Upwind_AdjFlow == UPWIND::ROE) cout << "Roe (with entropy fix = "<< EntropyFix_Coeff <<") solver for the adjoint inviscid terms."<< endl; if (MUSCL_AdjFlow) { PrintLimiterInfo(Kind_SlopeLimit_AdjFlow); } else { @@ -6687,7 +6689,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { if ((Kind_Solver == MAIN_SOLVER::ADJ_RANS) && (!Frozen_Visc_Cont)) { if (Kind_ConvNumScheme_AdjTurb == SPACE_UPWIND) { - if (Kind_Upwind_Turb == SCALAR_UPWIND) cout << "Scalar upwind solver for the adjoint turbulence model." << endl; + if (Kind_Upwind_Turb == UPWIND::SCALAR_UPWIND) cout << "Scalar upwind solver for the adjoint turbulence model." << endl; if (MUSCL_AdjTurb) { PrintLimiterInfo(Kind_SlopeLimit_AdjTurb); } else { @@ -6746,10 +6748,11 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "Discontinuous Galerkin Finite element solver" << endl; switch( Riemann_Solver_FEM ) { - case ROE: cout << "Roe (with entropy fix) solver for inviscid fluxes over the faces" << endl; break; - case LAX_FRIEDRICH: cout << "Lax-Friedrich solver for inviscid fluxes over the faces" << endl; break; - case AUSM: cout << "AUSM solver inviscid fluxes over the faces" << endl; break; - case HLLC: cout << "HLLC solver inviscid fluxes over the faces" << endl; break; + case UPWIND::ROE: cout << "Roe (with entropy fix) solver for inviscid fluxes over the faces" << endl; break; + case UPWIND::LAX_FRIEDRICH: cout << "Lax-Friedrich solver for inviscid fluxes over the faces" << endl; break; + case UPWIND::AUSM: cout << "AUSM solver inviscid fluxes over the faces" << endl; break; + case UPWIND::HLLC: cout << "HLLC solver inviscid fluxes over the faces" << endl; break; + default: break; } if(Kind_Solver != MAIN_SOLVER::FEM_EULER && Kind_Solver != MAIN_SOLVER::DISC_ADJ_FEM_EULER) { @@ -8287,7 +8290,7 @@ unsigned short CConfig::GetContainerPosition(unsigned short val_eqsystem) { } void CConfig::SetKind_ConvNumScheme(unsigned short val_kind_convnumscheme, - unsigned short val_kind_centered, unsigned short val_kind_upwind, + CENTERED val_kind_centered, UPWIND val_kind_upwind, LIMITER val_kind_slopelimit, bool val_muscl, unsigned short val_kind_fem) { Kind_ConvNumScheme = val_kind_convnumscheme; @@ -8354,7 +8357,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver, SetSpeciesParam(); if (val_system == RUNTIME_HEAT_SYS) { - SetKind_ConvNumScheme(Kind_ConvNumScheme_Heat, NONE, NONE, LIMITER::NONE, NONE, NONE); + SetKind_ConvNumScheme(Kind_ConvNumScheme_Heat, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, NONE); SetKind_TimeIntScheme(Kind_TimeIntScheme_Heat); } break; @@ -8370,7 +8373,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver, SetKind_TimeIntScheme(Kind_TimeIntScheme_Turb); } if (val_system == RUNTIME_HEAT_SYS) { - SetKind_ConvNumScheme(Kind_ConvNumScheme_Heat, NONE, NONE, LIMITER::NONE, NONE, NONE); + SetKind_ConvNumScheme(Kind_ConvNumScheme_Heat, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, NONE); SetKind_TimeIntScheme(Kind_TimeIntScheme_Heat); } break; @@ -8403,7 +8406,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver, break; case MAIN_SOLVER::HEAT_EQUATION: if (val_system == RUNTIME_HEAT_SYS) { - SetKind_ConvNumScheme(NONE, NONE, NONE, LIMITER::NONE, NONE, NONE); + SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, NONE); SetKind_TimeIntScheme(Kind_TimeIntScheme_Heat); } break; @@ -8413,7 +8416,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver, Current_DynTime = static_cast(TimeIter)*Delta_DynTime; if (val_system == RUNTIME_FEA_SYS) { - SetKind_ConvNumScheme(NONE, NONE, NONE, LIMITER::NONE , NONE, NONE); + SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE , NONE, NONE); SetKind_TimeIntScheme(NONE); } break; diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index abd45931ebc6..d493ac5c3659 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -55,6 +55,7 @@ class CNumerics { su2double Gas_Constant; /*!< \brief Gas constant. */ su2double Prandtl_Lam; /*!< \brief Laminar Prandtl's number. */ su2double Prandtl_Turb; /*!< \brief Turbulent Prandtl's number. */ + su2double MassFlux; /*!< \brief Mass flux across edge. */ su2double *Proj_Flux_Tensor; /*!< \brief Flux tensor projected in a direction. */ su2double **tau; /*!< \brief Viscous stress tensor. */ @@ -191,6 +192,8 @@ class CNumerics { bool nemo; /*!< \brief Flag for NEMO problems */ + bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ + public: /*! * \brief Return type used in some "ComputeResidual" overloads to give a @@ -1608,6 +1611,18 @@ class CNumerics { * \param[in] SolverSPvals - Struct holding the values. */ virtual void SetStreamwisePeriodicValues(const StreamwisePeriodicValues SolverSPvals) { } + + /*! + * \brief SetMassFlux + * \param[in] val_MassFlux: Mass flux across the edge + */ + inline void SetMassFlux(const su2double val_MassFlux) { MassFlux = val_MassFlux; } + + /*! + * \brief Obtain information on bounded scalar problem + * \return is_bounded_scalar : scalar solver uses bounded scalar convective transport + */ + inline bool GetBoundedScalar() const { return bounded_scalar;} }; /*! diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp index 0c06703588ff..0bee88104253 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp @@ -115,25 +115,32 @@ class CUpwScalar : public CNumerics { } AD::SetPreaccIn(&V_i[idx.Velocity()], nDim); AD::SetPreaccIn(&V_j[idx.Velocity()], nDim); + AD::SetPreaccIn(V_i[idx.Density()]); + AD::SetPreaccIn(V_j[idx.Density()]); + AD::SetPreaccIn(MassFlux); ExtraADPreaccIn(); - su2double q_ij = 0.0; - if (dynamic_grid) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - su2double Velocity_i = V_i[iDim + idx.Velocity()] - GridVel_i[iDim]; - su2double Velocity_j = V_j[iDim + idx.Velocity()] - GridVel_j[iDim]; - q_ij += 0.5 * (Velocity_i + Velocity_j) * Normal[iDim]; - } + if (bounded_scalar) { + a0 = fmax(0.0, MassFlux) / V_i[idx.Density()]; + a1 = fmin(0.0, MassFlux) / V_j[idx.Density()]; } else { - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - q_ij += 0.5 * (V_i[iDim + idx.Velocity()] + V_j[iDim + idx.Velocity()]) * Normal[iDim]; + su2double q_ij = 0.0; + if (dynamic_grid) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + su2double Velocity_i = V_i[iDim + idx.Velocity()] - GridVel_i[iDim]; + su2double Velocity_j = V_j[iDim + idx.Velocity()] - GridVel_j[iDim]; + q_ij += 0.5 * (Velocity_i + Velocity_j) * Normal[iDim]; + } + } else { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + q_ij += 0.5 * (V_i[iDim + idx.Velocity()] + V_j[iDim + idx.Velocity()]) * Normal[iDim]; + } } + a0 = fmax(0.0, q_ij); + a1 = fmin(0.0, q_ij); } - a0 = 0.5 * (q_ij + fabs(q_ij)); - a1 = 0.5 * (q_ij - fabs(q_ij)); - FinishResidualCalc(config); AD::SetPreaccOut(Flux, nVar); diff --git a/SU2_CFD/include/numerics/species/species_convection.hpp b/SU2_CFD/include/numerics/species/species_convection.hpp index 3308a8cd6916..64f296faff86 100644 --- a/SU2_CFD/include/numerics/species/species_convection.hpp +++ b/SU2_CFD/include/numerics/species/species_convection.hpp @@ -51,14 +51,12 @@ class CUpwSca_Species final : public CUpwScalar { using Base::ScalarVar_i; using Base::ScalarVar_j; using Base::idx; + using Base::bounded_scalar; /*! * \brief Adds any extra variables to AD */ - void ExtraADPreaccIn() override { - AD::SetPreaccIn(V_i[idx.Density()]); - AD::SetPreaccIn(V_j[idx.Density()]); - }; + void ExtraADPreaccIn() override {} /*! * \brief Species transport specific steps in the ComputeResidual method @@ -83,5 +81,5 @@ class CUpwSca_Species final : public CUpwScalar { * \param[in] config - Definition of the particular problem. */ CUpwSca_Species(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) - : CUpwScalar(val_nDim, val_nVar, config) {} + : CUpwScalar(val_nDim, val_nVar, config) { bounded_scalar = config->GetBounded_Species(); } }; diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp index b107cdff54a5..652a0c907b4d 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp @@ -48,6 +48,7 @@ class CUpwSca_TurbSA final : public CUpwScalar { using Base::ScalarVar_i; using Base::ScalarVar_j; using Base::implicit; + using Base::bounded_scalar; /*! * \brief Adds any extra variables to AD. @@ -75,7 +76,7 @@ class CUpwSca_TurbSA final : public CUpwScalar { * \param[in] config - Definition of the particular problem. */ CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) - : CUpwScalar(val_nDim, val_nVar, config) {} + : CUpwScalar(val_nDim, val_nVar, config) { bounded_scalar = config->GetBounded_Turb(); } }; /*! @@ -100,14 +101,12 @@ class CUpwSca_TurbSST final : public CUpwScalar { using Base::ScalarVar_j; using Base::implicit; using Base::idx; + using Base::bounded_scalar; /*! * \brief Adds any extra variables to AD */ - void ExtraADPreaccIn() override { - AD::SetPreaccIn(V_i[idx.Density()]); - AD::SetPreaccIn(V_j[idx.Density()]); - } + void ExtraADPreaccIn() override {} /*! * \brief SST specific steps in the ComputeResidual method @@ -116,7 +115,6 @@ class CUpwSca_TurbSST final : public CUpwScalar { void FinishResidualCalc(const CConfig* config) override { Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0]; Flux[1] = a0*V_i[idx.Density()]*ScalarVar_i[1] + a1*V_j[idx.Density()]*ScalarVar_j[1]; - if (implicit) { Jacobian_i[0][0] = a0; Jacobian_i[0][1] = 0.0; Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = a0; @@ -134,5 +132,5 @@ class CUpwSca_TurbSST final : public CUpwScalar { * \param[in] config - Definition of the particular problem. */ CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) - : CUpwScalar(val_nDim, val_nVar, config) {} + : CUpwScalar(val_nDim, val_nVar, config) { bounded_scalar = config->GetBounded_Turb(); } }; diff --git a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp index 805133920567..c967217c854d 100644 --- a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp +++ b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp @@ -41,9 +41,11 @@ template CNumericsSIMD* createUpwindIdealNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) { CNumericsSIMD* obj = nullptr; switch (config.GetKind_Upwind_Flow()) { - case ROE: + case UPWIND::ROE: obj = new CRoeScheme(config, iMesh, turbVars); break; + default: + break; } return obj; } @@ -62,19 +64,19 @@ CNumericsSIMD* createUpwindGeneralNumerics(const CConfig& config, int iMesh, con template CNumericsSIMD* createCenteredNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) { CNumericsSIMD* obj = nullptr; - switch ((iMesh==MESH_0)? config.GetKind_Centered_Flow() : LAX) { - case NO_CENTERED: + switch ((iMesh==MESH_0)? config.GetKind_Centered_Flow() : CENTERED::LAX) { + case CENTERED::NONE: break; - case LAX: + case CENTERED::LAX: obj = new CLaxScheme(config, iMesh, turbVars); break; - case JST: + case CENTERED::JST: obj = new CJSTScheme(config, iMesh, turbVars); break; - case JST_KE: + case CENTERED::JST_KE: obj = new CJSTkeScheme(config, iMesh, turbVars); break; - case JST_MAT: + case CENTERED::JST_MAT: obj = new CJSTmatScheme(config, iMesh, turbVars); break; } @@ -117,6 +119,8 @@ CNumericsSIMD* createNumerics(const CConfig& config, int iMesh, const CVariable* obj = createCenteredNumerics >(config, iMesh, turbVars); } break; + default: + break; } return obj; diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 4f005f1f4ad4..0b52f52d90b2 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -55,6 +55,8 @@ class CFVMFlowSolverBase : public CSolver { unsigned long omp_chunk_size; /*!< \brief Chunk size used in light point loops. */ + su2activevector EdgeMassFluxes; /*!< \brief Mass fluxes across each edge, for discretization of transported scalars. */ + /*! * \brief Utility to set the value of a member variables safely, and so that the new values are seen by all threads. * \param[in] lhsRhsPairs - Pairs of destination and source e.g. a,0,b,-1. @@ -2401,4 +2403,10 @@ class CFVMFlowSolverBase : public CSolver { inline su2double GetEddyViscWall(unsigned short val_marker, unsigned long val_vertex) const final { return EddyViscWall[val_marker][val_vertex]; } + + /*! + * \brief Get the mass fluxes across the edges (computed and stored during the discretization of convective fluxes). + */ + inline const su2activevector* GetEdgeMassFluxes() const final { return &EdgeMassFluxes; } + }; diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index e5e106112965..5301fef5ad3f 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -31,6 +31,8 @@ #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../variables/CScalarVariable.hpp" +#include "../variables/CFlowVariable.hpp" +#include "../variables/CPrimitiveIndices.hpp" #include "CSolver.hpp" /*! @@ -50,14 +52,15 @@ class CScalarSolver : public CSolver { unsigned long omp_chunk_size; /*!< \brief Chunk size used in light point loops. */ - su2double lowerlimit[MAXNVAR]; /*!< \brief contains lower limits for turbulence variables. Note that ::min() - returns the smallest positive value for floats. */ - su2double upperlimit[MAXNVAR]; /*!< \brief contains upper limits for turbulence variables. */ + su2double lowerlimit[MAXNVAR]; /*!< \brief contains lower limits for scalar variables. */ + su2double upperlimit[MAXNVAR]; /*!< \brief contains upper limits for scalar variables. */ su2double Solution_Inf[MAXNVAR]; /*!< \brief Far-field solution. */ const bool Conservative; /*!< \brief Transported Variable is conservative. Solution has to be multiplied with rho. */ + const CPrimitiveIndices prim_idx; /*!< \brief Indices of the primitive flow variables. */ + vector > SlidingState; // vector of matrix of pointers... inner dim alloc'd elsewhere (welcome, to the twilight zone) vector > SlidingStateNodes; @@ -86,7 +89,7 @@ class CScalarSolver : public CSolver { inline CVariable* GetBaseClassPointerToNodes() final { return nodes; } /*! - * \brief Compute the viscous flux for the turbulent equation at a particular edge. + * \brief Compute the viscous flux for the scalar equation at a particular edge. * \tparam SolverSpecificNumericsFunc - lambda-function, that implements solver specific contributions to numerics. * \note The functor has to implement (iPoint, jPoint) * \param[in] iEdge - Edge for which we want to compute the flux @@ -100,7 +103,7 @@ class CScalarSolver : public CSolver { CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, CConfig* config) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes(); + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); /*--- Points in edge ---*/ @@ -177,7 +180,7 @@ class CScalarSolver : public CSolver { su2double* PrimVar_i = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); - auto Jacobian_i = Jacobian.GetBlock(iPoint,iPoint); + auto* Jacobian_i = Jacobian.GetBlock(iPoint, iPoint); /*--- Loop over the nDonorVertexes and compute the averaged flux ---*/ @@ -208,6 +211,12 @@ class CScalarSolver : public CSolver { if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &PrimVar_j[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, true, density, velocity, Normal)); + } + auto residual = conv_numerics->ComputeResidual(config); /*--- Accumulate the residuals to compute the average ---*/ @@ -255,6 +264,20 @@ class CScalarSolver : public CSolver { } } + /*! + * \brief Applies a convective flux correction to negate the effects of flow divergence at a BC node. + * \note This function should be used for nodes that are part of a boundary marker, it computes a mass flux + * from density and velocity at the node, and the outward-pointing normal (-1 * normal of vertex). + * \return The mass flux. + */ + inline su2double BoundedScalarBCFlux(unsigned long iPoint, bool implicit, const su2double& density, + const su2double* velocity, const su2double* normal) { + const su2double edgeMassFlux = density * GeometryToolbox::DotProduct(nDim, velocity, normal); + LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -edgeMassFlux); + if (implicit) Jacobian.AddVal2Diag(iPoint, -edgeMassFlux); + return edgeMassFlux; + } + /*! * \brief Gradient and Limiter computation. * \param[in] geometry - Geometrical definition of the problem. @@ -265,7 +288,7 @@ class CScalarSolver : public CSolver { private: /*! - * \brief Compute the viscous flux for the turbulent equation at a particular edge. + * \brief Compute the viscous flux for the scalar equation at a particular edge. * \param[in] iEdge - Edge for which we want to compute the flux * \param[in] geometry - Geometrical definition of the problem. * \param[in] solver_container - Container vector with all the solutions. @@ -319,8 +342,20 @@ class CScalarSolver : public CSolver { * \param[in] config - Definition of the particular problem. * \param[in] iMesh - Index of the mesh in multigrid computations. */ - void Upwind_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config, - unsigned short iMesh) override; + void Upwind_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, + CConfig* config, unsigned short iMesh) override; + + /*! + * \brief Impose the Far Field 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_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) final; /*! * \brief Impose the Symmetry Plane boundary condition. diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 15752aa42a54..099e22680698 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -31,7 +31,9 @@ template CScalarSolver::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative) - : CSolver(), Conservative(conservative) { + : CSolver(), Conservative(conservative), + prim_idx(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE, + config->GetNEMOProblem(), geometry->GetnDim(), config->GetnSpecies()) { nMarker = config->GetnMarker_All(); /*--- Store the number of vertices on each marker for deallocation later ---*/ @@ -133,10 +135,14 @@ void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE) && (config->GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE); auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + const auto& edgeMassFluxes = *(solver_container[FLOW_SOL]->GetEdgeMassFluxes()); /*--- Pick one numerics object per thread. ---*/ auto* numerics = numerics_container[CONV_TERM + omp_get_thread_num() * MAX_TERMS]; + /*--- Apply scalar advection correction terms for bounded scalar problems ---*/ + const bool bounded_scalar = numerics->GetBoundedScalar(); + /*--- Static arrays of MUSCL-reconstructed flow primitives and turbulence variables (thread safety). ---*/ su2double solution_i[MAXNVAR] = {0.0}, flowPrimVar_i[MAXNVARFLOW] = {0.0}; su2double solution_j[MAXNVAR] = {0.0}, flowPrimVar_j[MAXNVARFLOW] = {0.0}; @@ -192,8 +198,10 @@ void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** Vector_ij[iDim] = 0.5 * (Coord_j[iDim] - Coord_i[iDim]); } - if (musclFlow) { - /*--- Reconstruct mean flow primitive variables. ---*/ + if (musclFlow && !bounded_scalar) { + /*--- Reconstruct mean flow primitive variables, note that in bounded scalar mode this is + * not necessary because the edge mass flux is read directly from the flow solver, instead + * of being computed from the primitive flow variables. ---*/ auto Gradient_i = flowNodes->GetGradient_Reconstruction(iPoint); auto Gradient_j = flowNodes->GetGradient_Reconstruction(jPoint); @@ -204,7 +212,8 @@ void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** } for (iVar = 0; iVar < solver_container[FLOW_SOL]->GetnPrimVarGrad(); iVar++) { - su2double Project_Grad_i = 0.0, Project_Grad_j = 0.0; + su2double Project_Grad_i = 0.0; + su2double Project_Grad_j = 0.0; for (iDim = 0; iDim < nDim; iDim++) { Project_Grad_i += Vector_ij[iDim] * Gradient_i[iVar][iDim]; Project_Grad_j -= Vector_ij[iDim] * Gradient_j[iVar][iDim]; @@ -249,6 +258,13 @@ void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** } } + /*--- Convective flux ---*/ + su2double EdgeMassFlux = 0.0; + if (bounded_scalar) { + EdgeMassFlux = edgeMassFluxes[iEdge]; + numerics->SetMassFlux(EdgeMassFlux); + } + /*--- Update convective residual value ---*/ auto residual = numerics->ComputeResidual(config); @@ -262,6 +278,21 @@ void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** if (implicit) Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); } + /*--- Apply convective flux correction to negate the effects of flow divergence in case of incompressible flow. + * Note that for the bounded scalar model, we explicitly put div(v)=0. + * If the ReducerStrategy is used, the corrections need to be applied in a loop over nodes + * to avoid race conditions in accessing nodes shared by edges handled by different threads. ---*/ + + if (bounded_scalar && !ReducerStrategy) { + LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -EdgeMassFlux); + LinSysRes.AddBlock(jPoint, nodes->GetSolution(jPoint), EdgeMassFlux); + + if (implicit) { + Jacobian.AddVal2Diag(iPoint, -EdgeMassFlux); + Jacobian.AddVal2Diag(jPoint, EdgeMassFlux); + } + } + /*--- Viscous contribution. ---*/ Viscous_Residual(iEdge, geometry, solver_container, @@ -277,6 +308,26 @@ void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** if (ReducerStrategy) { SumEdgeFluxes(geometry); if (implicit) Jacobian.SetDiagonalAsColumnSum(); + + /*--- Bounded scalar correction that cannot be applied in the edge loop when using the ReducerStrategy. ---*/ + if (bounded_scalar) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + const auto* solution = nodes->GetSolution(iPoint); + su2double divergence = 0; + + for (auto iEdge : geometry->nodes->GetEdges(iPoint)) { + const auto sign = (iPoint == geometry->edges->GetNode(iEdge,0)) ? 1 : -1; + const su2double EdgeMassFlux = sign * edgeMassFluxes[iEdge]; + divergence += EdgeMassFlux; + LinSysRes.AddBlock(iPoint, solution, -EdgeMassFlux); + } + if (implicit) { + Jacobian.AddVal2Diag(iPoint, -divergence); + } + } + END_SU2_OMP_FOR + } } } @@ -311,6 +362,67 @@ void CScalarSolver::BC_Periodic(CGeometry* geometry, CSolver** sol } } +template +void CScalarSolver::BC_Far_Field(CGeometry* geometry, CSolver** solver_container, + CNumerics* conv_numerics, CNumerics*, CConfig *config, + unsigned short val_marker) { + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + 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 infinity ---*/ + + auto V_infty = 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); + + /*--- Grid Movement ---*/ + + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + + conv_numerics->SetPrimitive(V_domain, V_infty); + + /*--- Set turbulent variable at the wall, and at infinity ---*/ + + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Solution_Inf); + + /*--- Set Normal (it is necessary to change the sign) ---*/ + + 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); + + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &V_infty[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal)); + } + + /*--- Compute residuals and Jacobians ---*/ + + auto residual = conv_numerics->ComputeResidual(config); + + /*--- Add residuals and Jacobians ---*/ + + LinSysRes.AddBlock(iPoint, residual); + if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } + } + END_SU2_OMP_FOR +} + template void CScalarSolver::PrepareImplicitIteration(CGeometry* geometry, CSolver** solver_container, CConfig* config) { diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index f5d268d40b5a..36c541117e54 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -30,6 +30,7 @@ #include "../../../Common/include/parallelization/mpi_structure.hpp" #include +#include #include #include #include @@ -2999,6 +3000,12 @@ class CSolver { */ inline virtual su2double GetEddyViscWall(unsigned short val_marker, unsigned long val_vertex) const { return 0; } + /*! + * \brief A virtual member + * \return The mass fluxes (from flow solvers) across the edges. + */ + inline virtual const su2activevector* GetEdgeMassFluxes() const { return nullptr; } + /*! * \brief A virtual member. * \return Value of the StrainMag_Max @@ -4323,7 +4330,7 @@ class CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. * \param[in] converged - Whether or not solution has converged. - */ + */ void SavelibROM(CGeometry *geometry, CConfig *config, bool converged); protected: diff --git a/SU2_CFD/include/solvers/CTransLMSolver.hpp b/SU2_CFD/include/solvers/CTransLMSolver.hpp index 36d29ee1a316..1ced793c9ac3 100644 --- a/SU2_CFD/include/solvers/CTransLMSolver.hpp +++ b/SU2_CFD/include/solvers/CTransLMSolver.hpp @@ -153,23 +153,7 @@ class CTransLMSolver final : public CTurbSolver { CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, - unsigned short val_marker) override; - - /*! - * \brief Impose the Far Field 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_Far_Field(CGeometry *geometry, - CSolver **solver_container, - CNumerics *conv_numerics, - CNumerics *visc_numerics, - CConfig *config, - unsigned short val_marker) override; + unsigned short val_marker) override; /*! * \brief Impose the inlet boundary condition. diff --git a/SU2_CFD/include/solvers/CTurbSASolver.hpp b/SU2_CFD/include/solvers/CTurbSASolver.hpp index 8868f8222635..8a09ad00cff9 100644 --- a/SU2_CFD/include/solvers/CTurbSASolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSASolver.hpp @@ -188,22 +188,6 @@ class CTurbSASolver final : public CTurbSolver { CConfig *config, unsigned short val_marker) override; - /*! - * \brief Impose the Far Field 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_Far_Field(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. diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp index cc881bfa6a31..c4b4625744d4 100644 --- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp @@ -169,22 +169,6 @@ class CTurbSSTSolver final : public CTurbSolver { CConfig *config, unsigned short val_marker) override; - /*! - * \brief Impose the Far Field 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_Far_Field(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. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 160dbbd11b81..4eed3037c645 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1259,7 +1259,7 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_Turb()) { - case NO_UPWIND: + case NO_CONVECTIVE: SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_TURB option.", CURRENT_FUNCTION); break; case SPACE_UPWIND : @@ -1350,7 +1350,7 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_Turb()) { - case NO_UPWIND: + case NONE: SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_TURB option.", CURRENT_FUNCTION); break; case SPACE_UPWIND : @@ -1682,8 +1682,8 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (incompressible) { /*--- Incompressible flow, use preconditioning method ---*/ switch (config->GetKind_Centered_Flow()) { - case LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLaxInc_Flow(nDim, nVar_Flow, config); break; - case JST : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentJSTInc_Flow(nDim, nVar_Flow, config); break; + case CENTERED::LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLaxInc_Flow(nDim, nVar_Flow, config); break; + case CENTERED::JST : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentJSTInc_Flow(nDim, nVar_Flow, config); break; default: SU2_MPI::Error("Invalid centered scheme or not implemented.\n Currently, only JST and LAX-FRIEDRICH are available for incompressible flows.", CURRENT_FUNCTION); break; @@ -1701,7 +1701,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (compressible) { /*--- Compressible flow ---*/ switch (config->GetKind_Upwind_Flow()) { - case ROE: + case UPWIND::ROE: if (ideal_gas) { for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { @@ -1717,62 +1717,62 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol } break; - case AUSM: + case UPWIND::AUSM: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwAUSM_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwAUSM_Flow(nDim, nVar_Flow, config); } break; - case AUSMPLUSUP: + case UPWIND::AUSMPLUSUP: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwAUSMPLUSUP_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwAUSMPLUSUP_Flow(nDim, nVar_Flow, config); } break; - case AUSMPLUSUP2: + case UPWIND::AUSMPLUSUP2: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwAUSMPLUSUP2_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwAUSMPLUSUP2_Flow(nDim, nVar_Flow, config); } break; - case TURKEL: + case UPWIND::TURKEL: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwTurkel_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwTurkel_Flow(nDim, nVar_Flow, config); } break; - case L2ROE: + case UPWIND::L2ROE: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwL2Roe_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwL2Roe_Flow(nDim, nVar_Flow, config); } break; - case LMROE: + case UPWIND::LMROE: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwLMRoe_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwLMRoe_Flow(nDim, nVar_Flow, config); } break; - case SLAU: + case UPWIND::SLAU: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwSLAU_Flow(nDim, nVar_Flow, config, roe_low_dissipation); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwSLAU_Flow(nDim, nVar_Flow, config, false); } break; - case SLAU2: + case UPWIND::SLAU2: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwSLAU2_Flow(nDim, nVar_Flow, config, roe_low_dissipation); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwSLAU2_Flow(nDim, nVar_Flow, config, false); } break; - case HLLC: + case UPWIND::HLLC: if (ideal_gas) { for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwHLLC_Flow(nDim, nVar_Flow, config); @@ -1787,14 +1787,14 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol } break; - case MSW: + case UPWIND::MSW: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwMSW_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwMSW_Flow(nDim, nVar_Flow, config); } break; - case CUSP: + case UPWIND::CUSP: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwCUSP_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwCUSP_Flow(nDim, nVar_Flow, config); @@ -1810,7 +1810,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (incompressible) { /*--- Incompressible flow, use artificial compressibility method ---*/ switch (config->GetKind_Upwind_Flow()) { - case FDS: + case UPWIND::FDS: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwFDSInc_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwFDSInc_Flow(nDim, nVar_Flow, config); @@ -1930,7 +1930,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (compressible) { /*--- Compressible flow ---*/ switch (config->GetKind_Centered_Flow()) { - case LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLax_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); break; + case CENTERED::LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLax_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); break; default: SU2_MPI::Error("Invalid centered scheme or not implemented.", CURRENT_FUNCTION); break; @@ -1948,35 +1948,35 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (compressible) { /*--- Compressible flow ---*/ switch (config->GetKind_Upwind_Flow()) { - case ROE: + case UPWIND::ROE: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwRoe_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwRoe_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); } break; - case AUSM: + case UPWIND::AUSM: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwAUSM_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwAUSM_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); } break; - case AUSMPLUSUP2: + case UPWIND::AUSMPLUSUP2: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwAUSMPLUSUP2_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwAUSMPLUSUP2_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); } break; - case MSW: + case UPWIND::MSW: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwMSW_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwMSW_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); } break; - case AUSMPWPLUS: + case UPWIND::AUSMPWPLUS: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwAUSMPWplus_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwAUSMPWplus_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); @@ -2020,41 +2020,41 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if ((fem_euler) || (fem_ns)) { switch (config->GetRiemann_Solver_FEM()) { - case ROE: - case LAX_FRIEDRICH: + case UPWIND::ROE: + case UPWIND::LAX_FRIEDRICH: /* Hard coded optimized implementation is used in the DG solver. No need to allocate the corresponding entry in numerics. */ break; - case AUSM: + case UPWIND::AUSM: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwAUSM_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwAUSM_Flow(nDim, nVar_Flow, config); } break; - case TURKEL: + case UPWIND::TURKEL: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwTurkel_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwTurkel_Flow(nDim, nVar_Flow, config); } break; - case HLLC: + case UPWIND::HLLC: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwHLLC_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwHLLC_Flow(nDim, nVar_Flow, config); } break; - case MSW: + case UPWIND::MSW: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwMSW_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwMSW_Flow(nDim, nVar_Flow, config); } break; - case CUSP: + case UPWIND::CUSP: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwCUSP_Flow(nDim, nVar_Flow, config); numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwCUSP_Flow(nDim, nVar_Flow, config); @@ -2170,8 +2170,8 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Compressible flow ---*/ switch (config->GetKind_Centered_AdjFlow()) { - case LAX : numerics[MESH_0][ADJFLOW_SOL][conv_term] = new CCentLax_AdjFlow(nDim, nVar_Adj_Flow, config); break; - case JST : numerics[MESH_0][ADJFLOW_SOL][conv_term] = new CCentJST_AdjFlow(nDim, nVar_Adj_Flow, config); break; + case CENTERED::LAX : numerics[MESH_0][ADJFLOW_SOL][conv_term] = new CCentLax_AdjFlow(nDim, nVar_Adj_Flow, config); break; + case CENTERED::JST : numerics[MESH_0][ADJFLOW_SOL][conv_term] = new CCentJST_AdjFlow(nDim, nVar_Adj_Flow, config); break; default: SU2_MPI::Error("Centered scheme not implemented.", CURRENT_FUNCTION); break; @@ -2193,7 +2193,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Compressible flow ---*/ switch (config->GetKind_Upwind_AdjFlow()) { - case ROE: + case UPWIND::ROE: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel][ADJFLOW_SOL][conv_term] = new CUpwRoe_AdjFlow(nDim, nVar_Adj_Flow, config); numerics[iMGlevel][ADJFLOW_SOL][conv_bound_term] = new CUpwRoe_AdjFlow(nDim, nVar_Adj_Flow, config); diff --git a/SU2_CFD/src/numerics/species/species_sources.cpp b/SU2_CFD/src/numerics/species/species_sources.cpp index 1198e6c757b1..50183ac85ea9 100644 --- a/SU2_CFD/src/numerics/species/species_sources.cpp +++ b/SU2_CFD/src/numerics/species/species_sources.cpp @@ -41,6 +41,7 @@ CSourceBase_Species::CSourceBase_Species(unsigned short val_nDim, unsigned short for (unsigned short iVar = 0; iVar < nVar; iVar++) { jacobian[iVar] = new su2double[nVar](); } + bounded_scalar = config->GetBounded_Species(); } CSourceBase_Species::~CSourceBase_Species() { @@ -67,13 +68,11 @@ CSourceAxisymmetric_Species::CSourceAxisymmetric_Species(unsigned short val_n template CNumerics::ResidualType<> CSourceAxisymmetric_Species::ComputeResidual(const CConfig* config) { - - /*--- Preaccumulation ---*/ AD::StartPreacc(); AD::SetPreaccIn(ScalarVar_i, nVar); - AD::SetPreaccIn(Volume); - + AD::SetPreaccIn(Volume); + if (incompressible) { AD::SetPreaccIn(V_i, nDim+6); } @@ -94,29 +93,31 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Species::ComputeResidual(const AD::SetPreaccIn(Coord_i[1]); AD::SetPreaccIn(Diffusion_Coeff_i, nVar); - AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); const su2double yinv = 1.0 / Coord_i[1]; const su2double Density_i = V_i[idx.Density()]; - + su2double Velocity_i[3]; for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_i[iDim] = V_i[idx.Velocity() + iDim]; - /*--- 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]; + /*--- Inviscid component of the source term. When div(v)=0, this term does not occur ---*/ + + if (!bounded_scalar) { + 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[iVar][iVar] -= yinv * Volume * Velocity_i[1]; + if (implicit) { + for (auto iVar = 0u; iVar < nVar; iVar++) { + jacobian[iVar][iVar] -= yinv * Volume * Velocity_i[1]; + } } } /*--- Add the viscous terms if necessary. ---*/ - + if (config->GetViscous()) { su2double Mass_Diffusivity_Tur = 0.0; if (turbulence) @@ -124,10 +125,11 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Species::ComputeResidual(const 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]; - } + } } + } - + AD::SetPreaccOut(residual, nVar); AD::EndPreacc(); diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp index 93bd29fa9bf8..1faf97175619 100644 --- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp @@ -951,7 +951,7 @@ void CAdjEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contai bool muscl = config->GetMUSCL_AdjFlow(); bool limiter = (config->GetKind_SlopeLimit_AdjFlow() != LIMITER::NONE); bool center = (config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED); - bool center_jst = (config->GetKind_Centered_AdjFlow() == JST); + bool center_jst = (config->GetKind_Centered_AdjFlow() == CENTERED::JST); bool fixed_cl = config->GetFixed_CL_Mode(); bool eval_dof_dcx = config->GetEval_dOF_dCX(); @@ -1034,7 +1034,7 @@ void CAdjEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co unsigned long iEdge, iPoint, jPoint; bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool jst_scheme = ((config->GetKind_Centered_AdjFlow() == JST) && (iMesh == MESH_0)); + bool jst_scheme = ((config->GetKind_Centered_AdjFlow() == CENTERED::JST) && (iMesh == MESH_0)); bool grid_movement = config->GetGrid_Movement(); for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index d2d34572ead9..a0a0c8f5b0c8 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -323,7 +323,7 @@ void CAdjNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); bool limiter = (config->GetKind_SlopeLimit_AdjFlow() != LIMITER::NONE); - bool center_jst = (config->GetKind_Centered_AdjFlow() == JST); + bool center_jst = (config->GetKind_Centered_AdjFlow() == CENTERED::JST); bool fixed_cl = config->GetFixed_CL_Mode(); bool eval_dof_dcx = config->GetEval_dOF_dCX(); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 2dac42dab922..ecf9cd26b885 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1536,17 +1536,17 @@ void CEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_con bool disc_adjoint = config->GetDiscrete_Adjoint(); bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); - bool center_jst = (config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0); - bool center_jst_ke = (config->GetKind_Centered_Flow() == JST_KE) && (iMesh == MESH_0); - bool center_jst_mat = (config->GetKind_Centered_Flow() == JST_MAT) && (iMesh == MESH_0); + bool center_jst = (config->GetKind_Centered_Flow() == CENTERED::JST) && (iMesh == MESH_0); + bool center_jst_ke = (config->GetKind_Centered_Flow() == CENTERED::JST_KE) && (iMesh == MESH_0); + bool center_jst_mat = (config->GetKind_Centered_Flow() == CENTERED::JST_MAT) && (iMesh == MESH_0); bool engine = ((config->GetnMarker_EngineInflow() != 0) || (config->GetnMarker_EngineExhaust() != 0)); bool actuator_disk = ((config->GetnMarker_ActDiskInlet() != 0) || (config->GetnMarker_ActDiskOutlet() != 0)); bool fixed_cl = config->GetFixed_CL_Mode(); unsigned short kind_row_dissipation = config->GetKind_RoeLowDiss(); bool roe_low_dissipation = (kind_row_dissipation != NO_ROELOWDISS) && - (config->GetKind_Upwind_Flow() == ROE || - config->GetKind_Upwind_Flow() == SLAU || - config->GetKind_Upwind_Flow() == SLAU2); + (config->GetKind_Upwind_Flow() == UPWIND::ROE || + config->GetKind_Upwind_Flow() == UPWIND::SLAU || + config->GetKind_Upwind_Flow() == UPWIND::SLAU2); /*--- Set the primitive variables ---*/ @@ -1739,14 +1739,14 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain const bool low_mach_corr = config->Low_Mach_Correction(); /*--- Use vectorization if the scheme supports it. ---*/ - if (config->GetKind_Upwind_Flow() == ROE && ideal_gas && !low_mach_corr) { + if (config->GetKind_Upwind_Flow() == UPWIND::ROE && ideal_gas && !low_mach_corr) { EdgeFluxResidual(geometry, solver_container, config); return; } const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - const bool roe_turkel = (config->GetKind_Upwind_Flow() == TURKEL); + const bool roe_turkel = (config->GetKind_Upwind_Flow() == UPWIND::TURKEL); const auto kind_dissipation = config->GetKind_RoeLowDiss(); const bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); @@ -2428,7 +2428,7 @@ void CEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, CCon return matrix; } - } precond(this, config->Low_Mach_Preconditioning() || (config->GetKind_Upwind_Flow() == TURKEL), nVar); + } precond(this, config->Low_Mach_Preconditioning() || (config->GetKind_Upwind_Flow() == UPWIND::TURKEL), nVar); PrepareImplicitIteration_impl(precond, geometry, config); } diff --git a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp index 2da6e045caf6..8779c86c3223 100644 --- a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp +++ b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp @@ -8923,7 +8923,7 @@ void CFEM_DG_EulerSolver::ComputeInviscidFluxesFace(CConfig *config /* Make a distinction between the several Riemann solvers. */ switch( config->GetRiemann_Solver_FEM() ) { - case ROE: { + case UPWIND::ROE: { /* Roe's approximate Riemann solver. Easier storage of the cut off value for the entropy correction. */ @@ -9168,7 +9168,7 @@ void CFEM_DG_EulerSolver::ComputeInviscidFluxesFace(CConfig *config /*------------------------------------------------------------------------*/ - case LAX_FRIEDRICH: { + case UPWIND::LAX_FRIEDRICH: { /* Local Lax-Friedrich (Rusanov) flux Make a distinction between two and three space dimensions in order to have the most efficient code. */ diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index b917e613670a..e7f238e7334c 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -222,6 +222,10 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned CommunicateInitialState(geometry, config); + /*--- Sizing edge mass flux array ---*/ + if (config->GetBounded_Scalar()) + EdgeMassFluxes.resize(geometry->GetnEdge()) = su2double(0.0); + /*--- Add the solver name (max 8 characters) ---*/ SolverName = "INC.FLOW"; @@ -900,7 +904,7 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_ const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); - const bool center_jst = (config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0); + const bool center_jst = (config->GetKind_Centered_Flow() == CENTERED::JST) && (iMesh == MESH_0); const bool outlet = (config->GetnMarker_Outlet() != 0); /*--- Set the primitive variables ---*/ @@ -1068,8 +1072,9 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co unsigned long iPoint, jPoint; - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool jst_scheme = ((config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0)); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool jst_scheme = ((config->GetKind_Centered_Flow() == CENTERED::JST) && (iMesh == MESH_0)); + const bool bounded_scalar = config->GetBounded_Scalar(); /*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of * variables, otherwise switch to the faster adjoint evaluation mode. ---*/ @@ -1117,6 +1122,8 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co auto residual = numerics->ComputeResidual(config); + if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0]; + /*--- Update residual value ---*/ if (ReducerStrategy) { @@ -1172,6 +1179,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont const bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); + const bool bounded_scalar = config->GetBounded_Scalar(); /*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of * variables, otherwise switch to the faster adjoint evaluation mode. ---*/ @@ -1290,6 +1298,8 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont auto residual = numerics->ComputeResidual(config); + if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0]; + /*--- Update residual value ---*/ if (ReducerStrategy) { @@ -1310,6 +1320,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont Viscous_Residual(iEdge, geometry, solver_container, numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); + } END_SU2_OMP_FOR } // end color loop @@ -2903,7 +2914,7 @@ void CIncEulerSolver::GetOutlet_Properties(CGeometry *geometry, CConfig *config, geometry->vertex[iMarker][iVertex]->GetNormal(Vector); - su2double AxiFactor = 1.0; + su2double AxiFactor = 1.0; if (axisymmetric) { if (geometry->nodes->GetCoord(iPoint, 1) > EPS) AxiFactor = 2.0*PI_NUMBER*geometry->nodes->GetCoord(iPoint, 1); @@ -2914,7 +2925,7 @@ void CIncEulerSolver::GetOutlet_Properties(CGeometry *geometry, CConfig *config, } } } - } + } Density = V_outlet[prim_idx.Density()]; @@ -3019,10 +3030,13 @@ void CIncEulerSolver::GetOutlet_Properties(CGeometry *geometry, CConfig *config, cout <<"Length (m): " << config->GetOutlet_Area(Outlet_TagBound) << "." << endl; } - cout << setprecision(5) << "Outlet Avg. Density (kg/m^3): " << config->GetOutlet_Density(Outlet_TagBound) * config->GetDensity_Ref() << endl; + cout << setprecision(5) << "Outlet Avg. Density (kg/m^3): " << config->GetOutlet_Density(Outlet_TagBound) * config->GetDensity_Ref() << endl; su2double Outlet_mDot = fabs(config->GetOutlet_MassFlow(Outlet_TagBound)) * config->GetDensity_Ref() * config->GetVelocity_Ref(); - cout << "Outlet mass flow (kg/s): "; cout << setprecision(5) << Outlet_mDot; - + su2double Outlet_mDot_Target = fabs(config->GetOutlet_Pressure(Outlet_TagBound)) / (config->GetDensity_Ref() * config->GetVelocity_Ref()); + cout << "Outlet mass flow (kg/s): " << setprecision(5) << Outlet_mDot << endl; + cout << "target mass flow (kg/s): " << setprecision(5) << Outlet_mDot_Target << endl; + su2double goal = 100.0*Outlet_mDot/Outlet_mDot_Target; + cout << "Target achieved:" << setprecision(5) << goal << " % "<< endl; } } diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 11a40c3a8390..79063af583bd 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -309,7 +309,7 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, c if (species_model != SPECIES_MODEL::NONE && solver_container[SPECIES_SOL] != nullptr) { scalar = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint); } - + /*--- Incompressible flow, primitive variables --- */ bool physical = static_cast(nodes)->SetPrimVar(iPoint,eddy_visc, turb_ke, GetFluidModel(), scalar); @@ -733,7 +733,7 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe const su2double U_Plus = VelTangMod / U_Tau; /*--- Y+ defined by White & Christoph ---*/ - + const su2double kUp = kappa * U_Plus; // incompressible adiabatic result diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index a347597a4a3e..7283e4259649 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -248,8 +248,8 @@ void CNEMOEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); - bool center_jst = (config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0); - bool center_jst_ke = (config->GetKind_Centered_Flow() == JST_KE) && (iMesh == MESH_0); + bool center_jst = (config->GetKind_Centered_Flow() == CENTERED::JST) && (iMesh == MESH_0); + bool center_jst_ke = (config->GetKind_Centered_Flow() == CENTERED::JST_KE) && (iMesh == MESH_0); /*--- Set the primitive variables ---*/ ompMasterAssignBarrier(ErrorCounter,0); @@ -595,7 +595,7 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con if (!chk_err_j) Gamma_j = ComputeConsistentExtrapolation(GetFluidModel(), nSpecies, Primitive_j, dPdU_j, dTdU_j, dTvedU_j, Eve_j, Cvve_j); /*--- Recompute Conserved variables if Roe or MSW scheme ---*/ - if ((config->GetKind_Upwind_Flow() == ROE) || (config->GetKind_Upwind_Flow() == MSW)){ + if ((config->GetKind_Upwind_Flow() == UPWIND::ROE) || (config->GetKind_Upwind_Flow() == UPWIND::MSW)){ if (!chk_err_i) RecomputeConservativeVector(Conserved_i, Primitive_i); if (!chk_err_j) RecomputeConservativeVector(Conserved_j, Primitive_j); } diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index 74f1efb4552f..743ccb73f89e 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -345,6 +345,8 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) { string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_STAT(OMP_MIN_SIZE) @@ -393,16 +395,23 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &V_inlet[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal)); + } + /*--- 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); // Unfinished viscous contribution removed before right after d8a0da9a00. Further testing required. + } } END_SU2_OMP_FOR @@ -486,6 +495,9 @@ void CSpeciesSolver::SetUniformInlet(const CConfig* config, unsigned short iMark void CSpeciesSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) { + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_STAT(OMP_MIN_SIZE) @@ -511,8 +523,8 @@ void CSpeciesSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container, Jacobian.DeleteValsRowi(total_index); } } else { // weak BC - /*--- Allocate the value at the outlet ---*/ + /*--- Allocate the value at the outlet ---*/ auto V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ @@ -524,8 +536,8 @@ void CSpeciesSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container, conv_numerics->SetPrimitive(V_domain, V_outlet); /*--- Set the species variables. Here we use a Neumann BC such - that the species variable is copied from the interior of the - domain to the outlet before computing the residual. ---*/ + that the species variable is copied from the interior of the + domain to the outlet before computing the residual. ---*/ conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint)); @@ -538,16 +550,21 @@ void CSpeciesSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container, if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - /*--- Compute the residual using an upwind scheme ---*/ + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &V_outlet[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal)); + } + /*--- 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); // Unfinished viscous contribution removed before right after d8a0da9a00. Further testing required. + } } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index 2f15bf8ba590..4aa62c917b71 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -203,18 +203,18 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai const su2double mu = flowNodes->GetLaminarViscosity(iPoint); const su2double dist = geometry->nodes->GetWall_Distance(iPoint); su2double VorticityMag = GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)); - su2double StrainMag =flowNodes->GetStrainMag(iPoint); + su2double StrainMag =flowNodes->GetStrainMag(iPoint); VorticityMag = max(VorticityMag, 1e-12); StrainMag = max(StrainMag, 1e-12); // safety against division by zero const su2double Intermittecy = nodes->GetSolution(iPoint,0); const su2double Re_t = nodes->GetSolution(iPoint,1); - const su2double Re_v = rho*dist*dist*StrainMag/mu; + const su2double Re_v = rho*dist*dist*StrainMag/mu; const su2double omega = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,1); const su2double k = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); su2double Corr_Rec = 0.0; if(Re_t <= 1870){ - Corr_Rec = Re_t - (396.035e-02 + (-120.656e-04)*Re_t + (868.230e-06)*pow(Re_t, 2.) + Corr_Rec = Re_t - (396.035e-02 + (-120.656e-04)*Re_t + (868.230e-06)*pow(Re_t, 2.) +( -696.506e-09)*pow(Re_t, 3.) + (174.105e-12)*pow(Re_t, 4.)); } else { Corr_Rec = Re_t - ( 593.11 + (Re_t - 1870.0) * 0.482); @@ -262,7 +262,7 @@ void CTransLMSolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, void CTransLMSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { - + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); @@ -279,7 +279,7 @@ void CTransLMSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta SU2_OMP_FOR_DYN(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { - + /*--- Conservative variables w/o reconstruction ---*/ @@ -295,7 +295,7 @@ void CTransLMSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta numerics->SetScalarVar(turbNodes->GetSolution(iPoint), nullptr); numerics->SetScalarVarGradient(turbNodes->GetGradient(iPoint), nullptr); - /*--- Transition variables w/o reconstruction, and its gradient ---*/ + /*--- Transition variables w/o reconstruction, and its gradient ---*/ numerics->SetTransVar(nodes->GetSolution(iPoint), nullptr); numerics->SetTransVarGradient(nodes->GetGradient(iPoint), nullptr); @@ -307,7 +307,7 @@ void CTransLMSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Set distance to the surface ---*/ numerics->SetDistance(geometry->nodes->GetWall_Distance(iPoint), 0.0); - + /*--- Set vorticity and strain rate magnitude ---*/ numerics->SetVorticity(flowNodes->GetVorticity(iPoint), nullptr); @@ -316,7 +316,7 @@ void CTransLMSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Set coordnate (for debugging) ---*/ numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr); - + /*--- Compute the source term ---*/ auto residual = numerics->ComputeResidual(config); @@ -400,62 +400,6 @@ void CTransLMSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_co } - - -void CTransLMSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) - for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - 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 infinity ---*/ - - auto V_infty = 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); - - conv_numerics->SetPrimitive(V_domain, V_infty); - - /*--- Set turbulent variable at the wall, and at infinity ---*/ - - conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Solution_Inf); - - /*--- Set Normal (it is necessary to change the sign) ---*/ - - 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); - - /*--- Grid Movement ---*/ - - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); - - /*--- Compute residuals and Jacobians ---*/ - - auto residual = conv_numerics->ComputeResidual(config); - - /*--- Add residuals and Jacobians ---*/ - - LinSysRes.AddBlock(iPoint, residual); - if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - } - } - END_SU2_OMP_FOR - -} - void CTransLMSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); @@ -524,15 +468,12 @@ void CTransLMSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C } -void CTransLMSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, - CConfig *config, unsigned short val_marker) { - -BC_Far_Field(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); - +void CTransLMSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + BC_Far_Field(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); } - -void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, +void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, bool val_update_geo) { const string restart_filename = config->GetFilename(config->GetSolution_FileName(), "", val_iter); @@ -561,7 +502,7 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi const bool energy = config->GetEnergy_Equation(); const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); - if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; + if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; /*--- Load data from the restart into correct containers. ---*/ @@ -652,5 +593,5 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi Restart_Data = nullptr; } END_SU2_OMP_SAFE_GLOBAL_ACCESS - + } \ No newline at end of file diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 5183536d44db..2913591606a9 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -464,61 +464,6 @@ void CTurbSASolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_con } -void CTurbSASolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) - for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - 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 infinity ---*/ - - auto V_infty = 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); - - /*--- Grid Movement ---*/ - - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - - conv_numerics->SetPrimitive(V_domain, V_infty); - - /*--- Set turbulent variable at the wall, and at infinity ---*/ - - conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Solution_Inf); - - /*--- Set Normal (it is necessary to change the sign) ---*/ - - 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); - - /*--- Compute residuals and Jacobians ---*/ - - auto residual = conv_numerics->ComputeResidual(config); - - /*--- Add residuals and Jacobians ---*/ - - LinSysRes.AddBlock(iPoint, residual); - if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - - } - } - END_SU2_OMP_FOR - -} - void CTurbSASolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { @@ -572,6 +517,12 @@ void CTurbSASolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CN conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &V_inlet[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal)); + } + /*--- Compute the residual using an upwind scheme ---*/ auto residual = conv_numerics->ComputeResidual(config); @@ -610,7 +561,6 @@ void CTurbSASolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CN } } END_SU2_OMP_FOR - } void CTurbSASolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, @@ -658,6 +608,12 @@ void CTurbSASolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, C conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &V_outlet[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal)); + } + /*--- Compute the residual using an upwind scheme ---*/ auto residual = conv_numerics->ComputeResidual(config); @@ -696,7 +652,6 @@ void CTurbSASolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, C } } END_SU2_OMP_FOR - } void CTurbSASolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 16ec987d8b18..7b20fa1f44fe 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -28,7 +28,6 @@ #include "../../include/solvers/CTurbSSTSolver.hpp" #include "../../include/variables/CTurbSSTVariable.hpp" #include "../../include/variables/CFlowVariable.hpp" -#include "../../include/variables/CPrimitiveIndices.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" @@ -519,69 +518,11 @@ void CTurbSSTSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_co } -void CTurbSSTSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) - for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - 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 infinity ---*/ - - auto V_infty = 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); - - conv_numerics->SetPrimitive(V_domain, V_infty); - - /*--- Set turbulent variable at the wall, and at infinity ---*/ - - conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Solution_Inf); - - /*--- Set Normal (it is necessary to change the sign) ---*/ - - 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); - - /*--- Grid Movement ---*/ - - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); - - /*--- Compute residuals and Jacobians ---*/ - - auto residual = conv_numerics->ComputeResidual(config); - - /*--- Add residuals and Jacobians ---*/ - - LinSysRes.AddBlock(iPoint, residual); - if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - } - } - END_SU2_OMP_FOR - -} - void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - const CPrimitiveIndices prim_idx(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE, - config->GetNEMOProblem(), nDim, config->GetnSpecies()); - /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_STAT(OMP_MIN_SIZE) @@ -658,6 +599,12 @@ void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &V_inlet[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal)); + } + /*--- Compute the residual using an upwind scheme ---*/ auto residual = conv_numerics->ComputeResidual(config); @@ -701,7 +648,6 @@ void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C } END_SU2_OMP_FOR - } void CTurbSSTSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, @@ -750,6 +696,12 @@ void CTurbSSTSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + if (conv_numerics->GetBoundedScalar()) { + const su2double* velocity = &V_outlet[prim_idx.Velocity()]; + const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal)); + } + /*--- Compute the residual using an upwind scheme ---*/ auto residual = conv_numerics->ComputeResidual(config); @@ -792,7 +744,6 @@ void CTurbSSTSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, } } END_SU2_OMP_FOR - } diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 3f5c9684a20d..75833a37f67e 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1342,6 +1342,15 @@ def main(): species2_primitiveVenturi_mixingmodel.new_output = True test_list.append(species2_primitiveVenturi_mixingmodel) + # 2 species (1 eq) primitive venturi mixing using mixing model and bounded scalar transport + species2_primitiveVenturi_mixingmodel_boundedscalar = TestCase('species2_primitiveVenturi_mixingmodel_boundedscalar') + species2_primitiveVenturi_mixingmodel_boundedscalar.cfg_dir = "species_transport/venturi_primitive_3species" + species2_primitiveVenturi_mixingmodel_boundedscalar.cfg_file = "species2_primitiveVenturi_mixingmodel_boundedscalar.cfg" + species2_primitiveVenturi_mixingmodel_boundedscalar.test_iter = 50 + species2_primitiveVenturi_mixingmodel_boundedscalar.test_vals = [-5.419758, -4.490843, -4.496264, -5.724852, -0.120393, -5.693152, 5.000000, -1.766881, 5.000000, -4.958352, 5.000000, -2.150266, 0.000300, 0.000300, 0.000000, 0.000000] + species2_primitiveVenturi_mixingmodel_boundedscalar.new_output = True + test_list.append(species2_primitiveVenturi_mixingmodel_boundedscalar) + # 2 species (1 eq) primitive venturi mixing using mixing model including viscosity and thermal conductivity species2_primitiveVenturi_mixingmodel_viscosity = TestCase('species2_primitiveVenturi_mixingmodel_viscosity') species2_primitiveVenturi_mixingmodel_viscosity.cfg_dir = "species_transport/venturi_primitive_3species" @@ -1384,6 +1393,15 @@ def main(): species2_primitiveVenturi.new_output = True test_list.append(species2_primitiveVenturi) + # 2 species (1 eq) primitive venturi mixing with bounded scalar transport + species_primitiveVenturi_boundedscalar = TestCase('species2_primitiveVenturi_bounded_scalar') + species_primitiveVenturi_boundedscalar.cfg_dir = "species_transport/venturi_primitive_3species" + species_primitiveVenturi_boundedscalar.cfg_file = "species2_primitiveVenturi_boundedscalar.cfg" + species_primitiveVenturi_boundedscalar.test_iter = 50 + species_primitiveVenturi_boundedscalar.test_vals = [-5.297585, -4.397797, -4.377086, -5.593131, -1.011782, -5.623540, 5.000000, -1.775123, 5.000000, -4.086339, 5.000000, -2.080187, 0.000424, 0.000424, 0.000000, 0.000000] + species_primitiveVenturi_boundedscalar.new_output = True + test_list.append(species_primitiveVenturi_boundedscalar) + # 3 species (2 eq) primitive venturi mixing with inlet files. # Note that the residuals are exactly the same as for the non-inlet case which should be the case for a fresh inlet file. species3_primitiveVenturi_inletFile = TestCase('species3_primitiveVenturi_inletFile') diff --git a/TestCases/species_transport/venturi_primitive_3species/DAspecies3_primitiveVenturi.cfg b/TestCases/species_transport/venturi_primitive_3species/DAspecies3_primitiveVenturi.cfg index d90e8d749b40..16e3ef50617e 100644 --- a/TestCases/species_transport/venturi_primitive_3species/DAspecies3_primitiveVenturi.cfg +++ b/TestCases/species_transport/venturi_primitive_3species/DAspecies3_primitiveVenturi.cfg @@ -53,7 +53,7 @@ INLET_FILENAME= inlet_venturi.dat INC_INLET_TYPE= VELOCITY_INLET VELOCITY_INLET MARKER_INLET= ( gas_inlet, 300, 1.0, 1.0, 0.0, 0.0,\ air_axial_inlet, 300, 1.0, 0.0, -1.0, 0.0 ) -SPECIES_USE_STRONG_BC= NO +SPECIES_USE_STRONG_BC= YES MARKER_INLET_SPECIES= (gas_inlet, 0.5, 0.5,\ air_axial_inlet, 0.6, 0.0 ) % @@ -128,7 +128,7 @@ CONV_FILENAME= history MARKER_ANALYZE= outlet MARKER_ANALYZE_AVERAGE= MASSFLUX % -OUTPUT_FILES= RESTART_ASCII, PARAVIEW_MULTIBLOCK +OUTPUT_FILES= RESTART_ASCII, PARAVIEW VOLUME_OUTPUT= RESIDUAL, PRIMITIVE OUTPUT_WRT_FREQ= 1000 % diff --git a/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_boundedscalar.cfg b/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_boundedscalar.cfg new file mode 100644 index 000000000000..3055bfd07196 --- /dev/null +++ b/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_boundedscalar.cfg @@ -0,0 +1,132 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Species mixing with 2 species, i.e. 1 transport equations % +% Author: T. Kattmann % +% Institution: Bosch Thermotechniek B.V. % +% Date: 2021/10/14 % +% File Version 7.4.0 "Blackbird" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_RANS +KIND_TURB_MODEL= SST +% +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_DENSITY_INIT= 1.1766 +% +INC_VELOCITY_INIT= ( 1.00, 0.0, 0.0 ) +% +INC_ENERGY_EQUATION= YES +INC_TEMPERATURE_INIT= 300.0 +% +INC_NONDIM= DIMENSIONAL +% +% -------------------- FLUID PROPERTIES ------------------------------------- % +% +FLUID_MODEL= CONSTANT_DENSITY +% +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY +THERMAL_CONDUCTIVITY_CONSTANT= 0.0357 +% +PRANDTL_LAM= 0.72 +TURBULENT_CONDUCTIVITY_MODEL= NONE +PRANDTL_TURB= 0.90 +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.716E-5 +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( wall, 0.0 ) +MARKER_SYM= ( axis ) +% +INC_INLET_TYPE= VELOCITY_INLET VELOCITY_INLET +MARKER_INLET= ( gas_inlet, 300, 1.0, 1.0, 0.0, 0.0,\ + air_axial_inlet, 300, 1.0, 0.0, -1.0, 0.0 ) +SPECIES_USE_STRONG_BC= NO +MARKER_INLET_SPECIES= (gas_inlet, 1.0,\ + air_axial_inlet, 0.6 ) +% +INC_OUTLET_TYPE= PRESSURE_OUTLET +MARKER_OUTLET= ( outlet, 0.0) +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +% +CFL_NUMBER= 60 +CFL_REDUCTION_SPECIES= 1.0 +CFL_REDUCTION_TURB= 1.0 +% +ITER= 51 +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-8 +LINEAR_SOLVER_ITER= 5 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW = NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +KIND_SCALAR_MODEL= SPECIES_TRANSPORT +DIFFUSIVITY_MODEL= CONSTANT_DIFFUSIVITY +DIFFUSIVITY_CONSTANT= 0.001 +% +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR +MUSCL_SPECIES= NO +SLOPE_LIMITER_SPECIES = NONE +% +TIME_DISCRE_SPECIES= EULER_IMPLICIT +% +SPECIES_INIT= 1.0 +SPECIES_CLIPPING= YES +SPECIES_CLIPPING_MIN= 0.0 +SPECIES_CLIPPING_MAX= 1.0 +% +% -------------------- TURBULENT TRANSPORT ---------------------------------------% +% +CONV_NUM_METHOD_TURB= BOUNDED_SCALAR +MUSCL_TURB= NO +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE, RMS_SPECIES +CONV_RESIDUAL_MINVAL= -18 +CONV_STARTITER= 10 +% +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= primitiveVenturi.su2 +% +SCREEN_OUTPUT= INNER_ITER WALL_TIME \ + RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_TKE RMS_DISSIPATION RMS_SPECIES_0 \ + LINSOL_ITER LINSOL_RESIDUAL \ + LINSOL_ITER_TURB LINSOL_RESIDUAL_TURB \ + LINSOL_ITER_SPECIES LINSOL_RESIDUAL_SPECIES SURFACE_SPECIES_VARIANCE +SCREEN_WRT_FREQ_INNER= 10 +% +HISTORY_OUTPUT= RMS_RES FLOW_COEFF LINSOL SPECIES_COEFF SPECIES_COEFF_SURF +MARKER_ANALYZE= outlet gas_inlet air_axial_inlet +% +OUTPUT_FILES= RESTART_ASCII, PARAVIEW_MULTIBLOCK +VOLUME_OUTPUT= RESIDUAL, PRIMITIVE +OUTPUT_WRT_FREQ= 1000 +% +RESTART_SOL= NO +SOLUTION_FILENAME= solution +RESTART_FILENAME= restart.dat +% +WRT_PERFORMANCE= YES diff --git a/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_mixingmodel_boundedscalar.cfg b/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_mixingmodel_boundedscalar.cfg new file mode 100644 index 000000000000..4c668a8ca0a0 --- /dev/null +++ b/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_mixingmodel_boundedscalar.cfg @@ -0,0 +1,142 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Species mixing with 2 species, i.e. 1 transport equations % +% Including mixture dependent density, viscosity, thermal % +% conductivity % +% Author: Cristopher Morales Ubal % +% Institution: Eindhoven University of Technology % +% Date: 2022/06/15 % +% File Version 7.4.0 "Blackbird" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_RANS +KIND_TURB_MODEL= SST +% +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= VARIABLE +INC_DENSITY_INIT= 1.1766 +% +INC_VELOCITY_INIT= ( 1.00, 0.0, 0.0 ) +% +INC_ENERGY_EQUATION= YES +INC_TEMPERATURE_INIT= 300.0 +% +INC_NONDIM= DIMENSIONAL +% +% -------------------- FLUID PROPERTIES ------------------------------------- % +% +FLUID_MODEL= FLUID_MIXTURE +% +MOLECULAR_WEIGHT= 28.96, 16.043 +% +SPECIFIC_HEAT_CP = 1009.39, 1009.39 +% +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY +THERMAL_CONDUCTIVITY_CONSTANT= 0.0357, 0.0357 +% +PRANDTL_LAM= 0.72, 0.72 +TURBULENT_CONDUCTIVITY_MODEL= NONE +PRANDTL_TURB= 0.90, 0.90 +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +% +MU_CONSTANT= 1.716E-5, 1.716E-5 +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( wall, 0.0 ) +MARKER_SYM= ( axis ) +% +SPECIFIED_INLET_PROFILE= NO +INLET_FILENAME= inlet_venturi.dat +INC_INLET_TYPE= VELOCITY_INLET VELOCITY_INLET +MARKER_INLET= ( gas_inlet, 300, 1.0, 1.0, 0.0, 0.0,\ + air_axial_inlet, 300, 1.0, 0.0, -1.0, 0.0 ) +SPECIES_USE_STRONG_BC= NO +MARKER_INLET_SPECIES= (gas_inlet, 1.0,\ + air_axial_inlet, 0.6) +% +INC_OUTLET_TYPE= PRESSURE_OUTLET +MARKER_OUTLET= ( outlet, 0.0) +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +% +CFL_NUMBER= 60 +CFL_REDUCTION_SPECIES= 1.0 +CFL_REDUCTION_TURB= 1.0 +% +ITER= 1000 +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-8 +LINEAR_SOLVER_ITER= 5 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW = NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +KIND_SCALAR_MODEL= SPECIES_TRANSPORT +DIFFUSIVITY_MODEL= CONSTANT_DIFFUSIVITY +DIFFUSIVITY_CONSTANT= 0.001 +% +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR +MUSCL_SPECIES= NO +SLOPE_LIMITER_SPECIES = NONE +% +TIME_DISCRE_SPECIES= EULER_IMPLICIT +% +SPECIES_INIT= 1.0 +SPECIES_CLIPPING= YES +SPECIES_CLIPPING_MIN= 0.0 +SPECIES_CLIPPING_MAX= 1.0 +% +% -------------------- TURBULENT TRANSPORT ---------------------------------------% +% +CONV_NUM_METHOD_TURB= BOUNDED_SCALAR +MUSCL_TURB= NO +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE, RMS_SPECIES +CONV_RESIDUAL_MINVAL= -18 +CONV_STARTITER= 10 +% +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= primitiveVenturi.su2 +% +SCREEN_OUTPUT= INNER_ITER WALL_TIME \ + RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_TKE RMS_DISSIPATION RMS_SPECIES_0 \ + LINSOL_ITER LINSOL_RESIDUAL \ + LINSOL_ITER_TURB LINSOL_RESIDUAL_TURB \ + LINSOL_ITER_SPECIES LINSOL_RESIDUAL_SPECIES SURFACE_SPECIES_VARIANCE +SCREEN_WRT_FREQ_INNER= 10 +% +HISTORY_OUTPUT= RMS_RES FLOW_COEFF LINSOL SPECIES_COEFF SPECIES_COEFF_SURF +MARKER_ANALYZE= outlet gas_inlet air_axial_inlet +% +OUTPUT_FILES= RESTART_ASCII, PARAVIEW_MULTIBLOCK +VOLUME_OUTPUT= RESIDUAL, PRIMITIVE +OUTPUT_WRT_FREQ= 1000 +% +RESTART_SOL= NO +SOLUTION_FILENAME= solution +RESTART_FILENAME= restart.dat +% +WRT_PERFORMANCE= YES diff --git a/TestCases/vandv/rans/30p30n/config.cfg b/TestCases/vandv/rans/30p30n/config.cfg index b3dbee99ca29..da6749b88b36 100644 --- a/TestCases/vandv/rans/30p30n/config.cfg +++ b/TestCases/vandv/rans/30p30n/config.cfg @@ -76,7 +76,7 @@ NEWTON_KRYLOV_DPARAM= ( 0.0, 1e-20, -3, 1e-5 ) % r0, tp, rf, e % % ------------------------ CONVERGENCE CRITERIA ------------------------- % % -ITER= 10000 +ITER= 20 CONV_RESIDUAL_MINVAL= -11.5 % % --------------------------- INPUT / OUTPUT ---------------------------- % diff --git a/config_template.cfg b/config_template.cfg index 57835eccf687..5b76a40e403d 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -717,7 +717,7 @@ MARKER_INLET_SPECIES= (inlet, 0.5, ..., inlet2, 0.6, ...) % Use strong inlet and outlet BC in the species solver SPECIES_USE_STRONG_BC= NO % -% Convective numerical method for species transport (SCALAR_UPWIND) +% Convective numerical method for species transport (SCALAR_UPWIND, BOUNDED_SCALAR) CONV_NUM_METHOD_SPECIES= SCALAR_UPWIND % % Monotonic Upwind Scheme for Conservation Laws (TVD) in the species equations. @@ -1275,7 +1275,7 @@ KIND_MATRIX_COLORING= GREEDY_COLORING % -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% % -% Convective numerical method (SCALAR_UPWIND) +% Convective numerical method (SCALAR_UPWIND, BOUNDED_SCALAR) CONV_NUM_METHOD_TURB= SCALAR_UPWIND % % Time discretization (EULER_IMPLICIT, EULER_EXPLICIT)