Consistent mass fluxes between bulk flow and scalar solvers and boundedness correction (#1697)

* working commit for incompressible
* Added correction factor to scalar convective residual to negate the effects of nonzero flow divergence on the scalar solver residual
* Changed SCALAR_ADVECTION upwind option to BOUNDED_SCALAR, added weak inlet boundary condition for BOUNDED_SCALAR option for the species solver
Co-authored-by: Evert <evertbun@xs4all.nl>
Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>

Author: 3 people

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<ENUM_CENTERED>(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,15 +4672,15 @@ 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.
    * \note This value is obtained from the config file, and it is constant
    *       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,15 +4724,15 @@ 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.
    * \note This value is obtained from the config file, and it is constant
    *       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,15 +4776,30 @@ 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.
    * \note This value is obtained from the config file, and it is constant
    *       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,15 +4815,15 @@ 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.
    * \note This value is obtained from the config file, and it is constant
    *       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.

Common/include/option_structure.hpp

@@ -809,67 +809,69 @@ static const MapType<std::string, ENUM_GUST_DIR> 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<std::string, ENUM_CENTERED> 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<std::string, CENTERED> 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<std::string, ENUM_UPWIND> 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<std::string, UPWIND> 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)
 };
 
 /*!

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<string>& 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;
   }
 };