[WIP] Pressure Based Solver

Common/include/CConfig.hpp

@@ -135,7 +135,8 @@ class CConfig {
   Hold_GridFixed,           /*!< \brief Flag hold fixed some part of the mesh during the deformation. */
   Axisymmetric,             /*!< \brief Flag for axisymmetric calculations */
   Enable_Cuda,              /*!< \brief Flag for switching GPU computing*/
-  Integrated_HeatFlux;      /*!< \brief Flag for heat flux BC whether it deals with integrated values.*/
+  Integrated_HeatFlux,      /*!< \brief Flag for heat flux BC whether it deals with integrated values.*/
+  Pressure_based;           /*!< \brief FLag to check if we are using a pressure-based system.*/
   su2double Buffet_k;       /*!< \brief Sharpness coefficient for buffet sensor.*/
   su2double Buffet_lambda;  /*!< \brief Offset parameter for buffet sensor.*/
   su2double Damp_Engine_Inflow;   /*!< \brief Damping factor for the engine inlet. */
@@ -538,6 +539,7 @@ class CConfig {
   Kind_TimeIntScheme_AdjTurb,   /*!< \brief Time integration for the adjoint turbulence model. */
   Kind_TimeIntScheme_Species,   /*!< \brief Time integration for the species model. */
   Kind_TimeIntScheme_Heat,      /*!< \brief Time integration for the wave equations. */
+  Kind_TimeIntScheme_Poisson,   /*!< \brief Time integration for the poisson (pressure correction) equation/ */
   Kind_TimeStep_Heat,           /*!< \brief Time stepping method for the (fvm) heat equation. */
   n_Datadriven_files;
 
@@ -587,6 +589,11 @@ class CConfig {
   Kind_Upwind_Heat,             /*!< \brief Upwind scheme for the heat transfer model. */
   Kind_Upwind_Template;         /*!< \brief Upwind scheme for the template model. */
 
+  ENUM_PBITER             /*< \brief Pressure-based solver for incompressible flows. */
+  Kind_PBIter;
+  ENUM_INCOMP_SYSTEM
+  Kind_Incomp_System;
+
   bool MUSCL,              /*!< \brief MUSCL scheme (for the runtime eq. system). */
   MUSCL_Flow,              /*!< \brief MUSCL scheme for the flow equations.*/
   MUSCL_Turb,              /*!< \brief MUSCL scheme for the turbulence equations.*/
@@ -640,6 +647,7 @@ class CConfig {
   su2double Deform_Linear_Solver_Error;          /*!< \brief Min error of the linear solver for the implicit formulation. */
   su2double Linear_Solver_Smoother_Relaxation;   /*!< \brief Relaxation factor for iterative linear smoothers. */
   unsigned long Linear_Solver_Iter;              /*!< \brief Max iterations of the linear solver for the implicit formulation. */
+  unsigned long Poisson_Linear_Solver_Iter;      /*!< \brief Max iterations of the linear solver for the poisson solver*/
   unsigned long Deform_Linear_Solver_Iter;       /*!< \brief Max iterations of the linear solver for the implicit formulation. */
   unsigned long Linear_Solver_Restart_Frequency; /*!< \brief Restart frequency of the linear solver for the implicit formulation. */
   unsigned long Linear_Solver_Restart_Deflation; /*!< \brief Number of vectors used for deflated restarts. */
@@ -649,6 +657,8 @@ class CConfig {
   su2double SemiSpan;                   /*!< \brief Wing Semi span. */
   su2double MSW_Alpha;                  /*!< \brief Coefficient for blending states in the MSW scheme. */
   su2double Roe_Kappa;                  /*!< \brief Relaxation of the Roe scheme. */
+  su2double RCFactor;                   /*!< \brief Relaxation for the Rhie-Chow interpolation contribution. */
+  su2double Relaxation_Factor_PBFlow;   /*!< \brief Relaxation coefficient of the flow corrections in the PB solver. */
   su2double Relaxation_Factor_Adjoint;  /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */
   su2double Relaxation_Factor_CHT;      /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */
   su2double EntropyFix_Coeff;           /*!< \brief Entropy fix coefficient. */
@@ -3992,6 +4002,30 @@ class CConfig {
    */
   ENUM_REGIME GetKind_Regime(void) const { return Kind_Regime; }
 
+  /*!
+   * \brief Kind of incompressible solver formulation.
+   * \return Kind of incompressible solver.
+   */
+  ENUM_INCOMP_SYSTEM GetKind_Incomp_System(void) const { return Kind_Incomp_System; }
+
+  /*!
+   * \brief Kind of iteration used for pressure based iterations.
+   * \return Kind of iteration used for pressure based iterations.
+   */
+  ENUM_PBITER GetKind_PBIter(void) const { return Kind_PBIter; }
+
+  /*!
+   * \brief Set the kind of incompressible solver formulation that is used.
+   * \param[in] val_system - the type of system to use.
+   */
+  void SetIncomp_System(unsigned short val_system);
+
+  /*!
+   * \brief Set the pressure based iteration method.
+   * \param[in] val_PBIter - The iteration method for the pressure based solver.
+   */
+  void SetPBIter(unsigned short val_PBIter);
+
   /*!
    * \brief Governing equations of the flow (it can be different from the run time equation).
    * \param[in] val_zone - Zone where the soler is applied.
@@ -4366,6 +4400,12 @@ class CConfig {
    */
   unsigned long GetLinear_Solver_Iter(void) const { return Linear_Solver_Iter; }
 
+  /*!
+   * \brief Get max number of iterations of the linear solver for the poisson equation.
+   * \return Max number of iterations of the linear solver for the poisson equation.
+   */
+  unsigned long GetPoisson_Linear_Solver_Iter(void) const { return Poisson_Linear_Solver_Iter; }
+
   /*!
    * \brief Get max number of iterations of the linear solver for the implicit formulation.
    * \return Max number of iterations of the linear solver for the implicit formulation.
@@ -4399,6 +4439,18 @@ class CConfig {
    * \return Relaxation factor.
    */
   su2double GetLinear_Solver_Smoother_Relaxation(void) const { return Linear_Solver_Smoother_Relaxation; }
+
+  /*!
+   * \brief Get the relaxation coefficient of the flow correction for PB solver.
+   * \return relaxation coefficient of the flow correction for PB solver
+   */
+  su2double GetRelaxation_Factor_PBFlow(void) const { return Relaxation_Factor_PBFlow; }
+
+  /*!
+   * \brief Get the relaxation coefficient for the Rhie-Chow interpolation in the PB solver.
+   * \return relaxation coefficient of the Rhie-Chow interpolation.
+   */
+  su2double GetRCFactor(void) const  { return RCFactor; }
 
   /*!
    * \brief Get the relaxation factor for solution updates of adjoint solvers.
@@ -4655,6 +4707,15 @@ class CConfig {
    */
   unsigned short GetKind_TimeIntScheme(void) const { return Kind_TimeNumScheme; }
 
+  /*!
+   * \brief Get the kind of integration scheme (explicit or implicit)
+   *        for the poisson/pressure correction equations.
+   * \note This value is obtained from the config file, and it is constant
+   *       during the computation.
+   * \return Kind of integration scheme for the poisson/pressure correction equations.
+   */
+  unsigned short GetKind_TimeIntScheme_Poisson(void) const { return Kind_TimeIntScheme_Poisson; }
+
   /*!
    * \brief Get the kind of convective numerical scheme.
    * \note This is the information that the code will use, the method will

Common/include/option_structure.hpp

@@ -261,6 +261,7 @@ enum class MAIN_SOLVER {
   FEM_RANS,                    /*!< \brief Definition of the finite element Reynolds-averaged Navier-Stokes' (RANS) solver. */
   FEM_LES,                     /*!< \brief Definition of the finite element Large Eddy Simulation Navier-Stokes' (LES) solver. */
   MULTIPHYSICS,
+  POISSON_EQUATION,            /*!< \brief Definition of the Poisson equation solver. */
   NEMO_EULER,                  /*!< \brief Definition of the NEMO Euler solver. */
   NEMO_NAVIER_STOKES,          /*!< \brief Definition of the NEMO NS solver. */
 };
@@ -336,6 +337,33 @@ static const MapType<std::string, STRUCT_COMPRESS> MatComp_Map = {
   MakePair("NEARLY_INCOMPRESSIBLE", STRUCT_COMPRESS::NEARLY_INCOMP)
 };
 
+/*!
+ * \brief Type of incompressible solver
+ */
+enum class ENUM_INCOMP_SYSTEM {
+  DENSITY_BASED,		/*!< \brief Density-based. */
+  PRESSURE_BASED,       /*!< \brief Pressure-based. */
+};
+static const MapType<std::string, ENUM_INCOMP_SYSTEM> Incomp_Map = {
+ MakePair("DENSITY_BASED", ENUM_INCOMP_SYSTEM::DENSITY_BASED)
+ MakePair("PRESSURE_BASED", ENUM_INCOMP_SYSTEM::PRESSURE_BASED)
+};
+
+/*!
+ * \brief Type of iteration
+ */
+enum class ENUM_PBITER {
+  SIMPLE,		/*!< \brief SIMPLE algorithm. */
+  SIMPLEC,		/*!< \brief SIMPLE algorithm. */
+  PISO,         /*!< \brief PISO algorithm. */
+};
+
+static const MapType<std::string, ENUM_PBITER> PBIter_Map = {
+ MakePair("SIMPLE", ENUM_PBITER::SIMPLE)
+ MakePair("SIMPLEC", ENUM_PBITER::SIMPLEC)
+ MakePair("PISO", ENUM_PBITER::PISO)
+};
+
 /*!
  * \brief Types of interpolators
  */
@@ -476,6 +504,7 @@ enum RUNTIME_TYPE {
   RUNTIME_ADJRAD_SYS = 24,    /*!< \brief One-physics case, the code is solving the adjoint radiation model. */
   RUNTIME_SPECIES_SYS = 25,   /*!< \brief One-physics case, the code is solving the species model. */
   RUNTIME_ADJSPECIES_SYS = 26,/*!< \brief One-physics case, the code is solving the adjoint species model. */
+  RUNTIME_POISSON_SYS = 27,   /*!< \brief One-physics case, the code is solving the poisson equation. */
 };
 
  enum SOLVER_TYPE : const int {
@@ -494,6 +523,7 @@ enum RUNTIME_TYPE {
    ADJSPECIES_SOL=12, /*!< \brief Position of the adjoint of the species solver. */
    FEA_SOL=0,        /*!< \brief Position of the Finite Element flow solution in the solver container array. */
    ADJFEA_SOL=1,     /*!< \brief Position of the continuous adjoint Finite Element flow solution in the solver container array. */
+   POISSON_SOL=13,   /*!< \brief Position of the poisson solution in the solver container array */
    TEMPLATE_SOL=0,   /*!< \brief Position of the template solution. */
  };
 
@@ -825,7 +855,8 @@ enum class CENTERED {
   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. */
-  LD2             /*!< \brief Low-Dissipation Low-Dispersion (LD2) centered scheme. */
+  LD2,            /*!< \brief Low-Dissipation Low-Dispersion (LD2) centered scheme. */
+  CDS             /*!< \brief Central Difference Scheme used for pressure based solver. */
 };
 static const MapType<std::string, CENTERED> Centered_Map = {
   MakePair("NONE", CENTERED::NONE)
@@ -834,6 +865,7 @@ static const MapType<std::string, CENTERED> Centered_Map = {
   MakePair("JST_MAT", CENTERED::JST_MAT)
   MakePair("LAX-FRIEDRICH", CENTERED::LAX)
   MakePair("LD2", CENTERED::LD2)
+  MakePair("CDS", CENTERED::CDS)
 };
 
 
@@ -860,7 +892,8 @@ enum class UPWIND {
   AUSMPLUSUP,             /*!< \brief AUSM+ -up numerical method (All Speed) */
   AUSMPLUSUP2,            /*!< \brief AUSM+ -up2 numerical method (All Speed) */
   AUSMPLUSM,              /*!< \breif AUSM+M numerical method. (NEMO Only)*/
-  BOUNDED_SCALAR          /*!< \brief Scalar advection numerical method. */
+  BOUNDED_SCALAR,         /*!< \brief Scalar advection numerical method. */
+  UDS                     /*!< \brief Upwind Difference Scheme used for pressure based solver. */
 };
 static const MapType<std::string, UPWIND> Upwind_Map = {
   MakePair("NONE", UPWIND::NONE)
@@ -882,6 +915,7 @@ static const MapType<std::string, UPWIND> Upwind_Map = {
   MakePair("SLAU2", UPWIND::SLAU2)
   MakePair("FDS", UPWIND::FDS)
   MakePair("LAX-FRIEDRICH", UPWIND::LAX_FRIEDRICH)
+  MakePair("UDS", UPWIND::UDS)
 };
 
 /*!
@@ -2694,6 +2728,7 @@ enum PERIODIC_QUANTITIES {
   PERIODIC_LIM_PRIM_1 ,  /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */
   PERIODIC_LIM_PRIM_2 ,  /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */
   PERIODIC_IMPLICIT   ,  /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */
+  PERIODIC_PRESSURE   ,  /*!< \brief Corrected pressure communication. */
 };
 
 /*!
@@ -2727,6 +2762,9 @@ enum class MPI_QUANTITIES {
   MESH_DISPLACEMENTS   ,  /*!< \brief Mesh displacements at the interface. */
   SOLUTION_TIME_N      ,  /*!< \brief Solution at time n. */
   SOLUTION_TIME_N1     ,  /*!< \brief Solution at time n-1. */
+  PRESSURE_VAR         ,  /*!< \brief Primitive variable (pressure) communication. */
+  MASS_FUX             ,  /*!< \brief Primitive variable (mass flux) communication. */
+  MOM_COEFF            ,  
 };
 
 /*!

Common/src/CConfig.cpp

@@ -1137,6 +1137,11 @@ void CConfig::SetConfig_Options() {
   /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SA turbulence model options/corrections. \n Options: see \link SA_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumListOption("SA_OPTIONS", nSA_Options, SA_Options, SA_Options_Map);
 
+  /*!\brief KIND_INCOMP_SYSTEM \n DESCRIPTION: Incomp type \n OPTIONS: see \link Incomp_Map \endlink DEFAULT: NONE \ingroup Config*/
+  addEnumOption("KIND_INCOMP_SYSTEM", Kind_Incomp_System, Incomp_Map, ENUM_INCOMP_SYSTEM::DENSITY_BASED);
+  /*!\brief KIND_PB_ITER \n  DESCRIPTION: Kind_PBIter \n OPTIONS: see \link PBIter_Map \endlink \ingroup Config*/
+  addEnumOption("KIND_PB_ITER", Kind_PBIter, PBIter_Map, ENUM_PBITER::SIMPLE);
+
   /*!\brief ROUGHSST_OPTIONS \n DESCRIPTION: Specify type of boundary condition for rough walls for SST turbulence model. \n Options: see \link ROUGHSST_Options_Map \endlink \n DEFAULT: wilcox1998 \ingroup Config*/
   addEnumOption("KIND_ROUGHSST_MODEL", Kind_RoughSST_Model, RoughSST_Model_Map, ROUGHSST_MODEL::WILCOX1998);
   /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/
@@ -1880,6 +1885,8 @@ void CConfig::SetConfig_Options() {
   addEnumOption("TIME_DISCRE_HEAT", Kind_TimeIntScheme_Heat, Time_Int_Map, EULER_IMPLICIT);
   /* DESCRIPTION: Time discretization */
   addEnumOption("TIMESTEP_HEAT", Kind_TimeStep_Heat, Heat_TimeStep_Map, MINIMUM);
+  /* DESCRIPTION: Time discretization */
+  addEnumOption("TIME_DISCRE_POISSON", Kind_TimeIntScheme_Poisson, Time_Int_Map, EULER_IMPLICIT);
 
   /*!\par CONFIG_CATEGORY: Linear solver definition \ingroup Config*/
   /*--- Options related to the linear solvers ---*/
@@ -1894,6 +1901,8 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("LINEAR_SOLVER_ERROR", Linear_Solver_Error, 1E-6);
   /* DESCRIPTION: Maximum number of iterations of the linear solver for the implicit formulation */
   addUnsignedLongOption("LINEAR_SOLVER_ITER", Linear_Solver_Iter, 10);
+  /* DESCRIPTION: Maximum number of iterations of the poisson linear solver for the implicit formulation */
+  addUnsignedLongOption("POISSON_LINEAR_SOLVER_ITER", Poisson_Linear_Solver_Iter, 10);
   /* DESCRIPTION: Fill in level for the ILU preconditioner */
   addUnsignedShortOption("LINEAR_SOLVER_ILU_FILL_IN", Linear_Solver_ILU_n, 0);
   /* DESCRIPTION: Use level scheduling for OMP parallelization of the ILU preconditioner */
@@ -1908,6 +1917,10 @@ void CConfig::SetConfig_Options() {
   addUnsignedLongOption("LINEAR_SOLVER_PREC_THREADS", Linear_Solver_Prec_Threads, 0);
   /* DESCRIPTION: Use an inner linear solver. */
   addEnumOption("LINEAR_SOLVER_INNER", Kind_Linear_Solver_Inner, Inner_Linear_Solver_Map, LINEAR_SOLVER_INNER::NONE);
+  /* DESCRIPTION: Relaxation of the pressure based flow corrections */
+  addDoubleOption("RELAXATION_FACTOR_PBFLOW", Relaxation_Factor_PBFlow, 0.5);
+  /* DESCRIPTION: Relaxation of the Rhie Chow interpolation contribution in pressure based flow. */
+  addDoubleOption("RELAXATION_FACTOR_RHIECHOW", RCFactor, 0.0);
   /* DESCRIPTION: Relaxation factor for updates of adjoint variables. */
   addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0);
   /* DESCRIPTION: Relaxation of the CHT coupling */
@@ -8805,6 +8818,7 @@ unsigned short CConfig::GetContainerPosition(unsigned short val_eqsystem) {
     case RUNTIME_ADJSPECIES_SYS:return ADJSPECIES_SOL;
     case RUNTIME_ADJFEA_SYS:    return ADJFEA_SOL;
     case RUNTIME_RADIATION_SYS: return RAD_SOL;
+    case RUNTIME_POISSON_SYS:   return POISSON_SOL;
     case RUNTIME_MULTIGRID_SYS: return 0;
   }
   return 0;
@@ -8943,6 +8957,13 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
       }
       break;
 
+    case MAIN_SOLVER::POISSON_EQUATION:
+    if (val_system == RUNTIME_POISSON_SYS) {
+      SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, 0.0, NONE);
+      SetKind_TimeIntScheme(Kind_TimeIntScheme_Poisson); 
+    }
+    break;
+
     case MAIN_SOLVER::FEM_ELASTICITY:
     case MAIN_SOLVER::DISC_ADJ_FEM:
       if (val_system == RUNTIME_FEA_SYS) {