[WIP] Pressure Based Flow Solver V8

Common/include/CConfig.hpp

@@ -133,7 +133,8 @@ class CConfig {
   Sens_Remove_Sharp,        /*!< \brief Flag for removing or not the sharp edges from the sensitivity computation. */
   Hold_GridFixed,           /*!< \brief Flag hold fixed some part of the mesh during the deformation. */
   Axisymmetric,             /*!< \brief Flag for axisymmetric calculations */
-  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. */
@@ -528,6 +529,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;
   ENUM_DATADRIVEN_METHOD Kind_DataDriven_Method;       /*!< \brief Method used for datset regression in data-driven fluid models. */
@@ -578,6 +580,12 @@ 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 .*/
   MUSCL_Flow,              /*!< \brief MUSCL scheme for the flow equations.*/
   MUSCL_Turb,              /*!< \brief MUSCL scheme for the turbulence equations.*/
@@ -621,12 +629,15 @@ 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 implicit formulation. */
   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_Prec_Threads;      /*!< \brief Number of threads per rank for ILU and LU_SGS preconditioners. */
   unsigned short Linear_Solver_ILU_n;            /*!< \brief ILU fill=in level. */
   su2double SemiSpan;                   /*!< \brief Wing Semi span. */
   su2double Roe_Kappa;                  /*!< \brief Relaxation of the Roe scheme. */
+  su2double RCFactor;                   /*!< \brief Relaxation for Rhie Chow interpolation contribution. */
+  su2double Relaxation_Factor_PBFlow;   /*!< \brief Relaxation coefficient of the flow corrections in 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. */
@@ -864,6 +875,11 @@ class CConfig {
   su2double Schmidt_Number_Laminar;   /*!< \brief Laminar Schmidt number for mass diffusion.  */
   su2double Schmidt_Number_Turbulent; /*!< \brief Turbulent Schmidt number for mass diffusion.  */
   su2double *Constant_Lewis_Number;   /*!< \brief Different Lewis number for mass diffusion.  */
+
+  su2double Mu_ConstantND,         /*!< \brief Non-dimensional constant viscosity for ConstantViscosity model.  */
+  Kt_ConstantND,         /*!< \brief Non-dimensional constant thermal conductivity for ConstantConductivity model.  */
+  Mu_SND,                /*!< \brief Non-dimensional reference S for Sutherland model.  */
+  Mu_Temperature_RefND;  /*!< \brief Non-dimensional reference temperature for Sutherland model.  */
   array<su2double, N_POLY_COEFFS> CpPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for specific heat Cp. */
   array<su2double, N_POLY_COEFFS> MuPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */
   array<su2double, N_POLY_COEFFS> KtPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */
@@ -3867,6 +3883,33 @@ class CConfig {
    */
   ENUM_REGIME GetKind_Regime(void) const { return Kind_Regime; }
 
+  /*!
+   * \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.
+   * \return Governing equation that we are solving.
+   */
+  ENUM_INCOMP_SYSTEM GetKind_Incomp_System(void) const { return Kind_Incomp_System; }
+
+  /*!
+   * \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.
+   * \return Governing equation that we are solving.
+   */
+  ENUM_PBITER GetKind_PBIter(void) const { return Kind_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.
+   * \return Governing equation that we are solving.
+   */
+  void SetIncomp_System(unsigned short val_system);
+
+  /*!
+   * \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.
+   * \return Governing equation that we are solving.
+   */
+  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.
@@ -4173,6 +4216,31 @@ class CConfig {
    */
   const su2double* GetKt_PolyCoeffND(void) const { return KtPolyCoefficientsND.data(); }
 
+  /*!
+   * \brief Set the value of the non-dimensional constant viscosity.
+   */
+  void SetMu_ConstantND(su2double mu_const) { Mu_ConstantND = mu_const; }
+
+  /*!
+   * \brief Set the value of the non-dimensional thermal conductivity.
+   */
+  void SetKt_ConstantND(su2double kt_const) { Kt_ConstantND = kt_const; }
+
+  /*!
+   * \brief Set the value of the non-dimensional reference viscosity for Sutherland model.
+   */
+  // void SetMu_RefND(su2double mu_ref) { Mu_RefND = mu_ref; }
+
+  /*!
+   * \brief Set the value of the non-dimensional reference temperature for Sutherland model.
+   */
+  void SetMu_Temperature_RefND(su2double mu_Tref) { Mu_Temperature_RefND = mu_Tref; }
+
+  /*!
+   * \brief Set the value of the non-dimensional S for Sutherland model.
+   */
+  void SetMu_SND(su2double mu_s) { Mu_SND = mu_s; }
+
   /*!
    * \brief Set the temperature polynomial coefficient for specific heat Cp.
    * \param[in] val_coeff - Temperature polynomial coefficient for specific heat Cp.
@@ -4260,6 +4328,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 implicit formulation.
+   * \return Max number of iterations of the linear solver for the implicit formulation.
+   */
+  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.
@@ -4283,6 +4357,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 of the flow correction for PB solver.
+   * \return relaxation coefficient of the flow correction for PB solver
+   */
+  su2double GetRCFactor(void) const  { return RCFactor; }
 
   /*!
    * \brief Get the relaxation factor for solution updates of adjoint solvers.
@@ -4640,6 +4726,15 @@ class CConfig {
    */
   unsigned short GetKind_TimeIntScheme_Radiation(void) const { return Kind_TimeIntScheme_Radiation; }
 
+  /*!
+   * \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 integration scheme (explicit or implicit)
    *        for the template equations.

Common/include/geometry/CGeometry.hpp

@@ -765,6 +765,11 @@ class CGeometry {
   inline virtual void SetAvgTurboValue(CConfig* config, unsigned short val_iZone, unsigned short marker_flag,
                                        bool allocate) {}
 
+  /*!
+   * \brief Sets CG coordinates.
+   */
+  inline virtual void SetCoord_CG(void) {}
+
   /*!
    * \brief A virtual member.
    * \param[in] config - Definition of the particular problem.

Common/include/geometry/CPhysicalGeometry.hpp

@@ -416,6 +416,11 @@ class CPhysicalGeometry final : public CGeometry {
    */
   void SetAvgTurboGeoValues(const CConfig* donor_config, CGeometry* donor_geometry, unsigned short donorZone) override;
 
+  /*!
+   * \brief Set the center of gravity of the face, elements and edges.
+   */
+  void SetCoord_CG(void) override;
+
   /*!
    * \brief Set the edge structure of the control volume.
    * \param[in] config - Definition of the particular problem.

Common/include/option_structure.hpp

@@ -256,6 +256,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,
   NEMO_EULER,                  /*!< \brief Definition of the NEMO Euler solver. */
   NEMO_NAVIER_STOKES,          /*!< \brief Definition of the NEMO NS solver. */
 };
@@ -331,6 +332,36 @@ static const MapType<std::string, STRUCT_COMPRESS> MatComp_Map = {
   MakePair("NEARLY_INCOMPRESSIBLE", STRUCT_COMPRESS::NEARLY_INCOMP)
 };
 
+/*------New option for pressure-based system-------*/
+/*!
+ * \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)
+};
+
+/*------ End new option for pressure-based system-------*/
+
 /*!
  * \brief Types of interpolators
  */
@@ -470,6 +501,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 adjoint species model. */
 };
 
  enum SOLVER_TYPE : const int {
@@ -488,6 +520,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=12,   /*!< \brief Position of the template solution. */
    TEMPLATE_SOL=0,   /*!< \brief Position of the template solution. */
  };
 
@@ -818,14 +851,16 @@ enum class CENTERED {
   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. */
+  JST_KE,          /*!< \brief Kinetic Energy preserving Jameson-Smith-Turkel centered numerical method. */
+  CDS               /*!< \brief Used for pressure based solver. */
 };
 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)
+  MakePair("CDS", CENTERED::CDS)
 };
 
 
@@ -852,7 +887,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
 };
 static const MapType<std::string, UPWIND> Upwind_Map = {
   MakePair("NONE", UPWIND::NONE)
@@ -873,6 +909,7 @@ static const MapType<std::string, UPWIND> Upwind_Map = {
   MakePair("LMROE", UPWIND::LMROE)
   MakePair("SLAU2", UPWIND::SLAU2)
   MakePair("FDS", UPWIND::FDS)
+  MakePair("UDS", UPWIND::UDS)
   MakePair("LAX-FRIEDRICH", UPWIND::LAX_FRIEDRICH)
 };
 
@@ -1868,10 +1905,12 @@ static const MapType<std::string, INLET_TYPE> Inlet_Map = {
 enum class INC_OUTLET_TYPE {
   PRESSURE_OUTLET,    /*!< \brief Gauge pressure outlet for incompressible flow */
   MASS_FLOW_OUTLET,   /*!< \brief Mass flow outlet for incompressible flow. */
+  OPEN,                /*!< \brief Inlet-Outlet type. */
 };
 static const MapType<std::string, INC_OUTLET_TYPE> Inc_Outlet_Map = {
   MakePair("PRESSURE_OUTLET",  INC_OUTLET_TYPE::PRESSURE_OUTLET)
   MakePair("MASS_FLOW_OUTLET", INC_OUTLET_TYPE::MASS_FLOW_OUTLET)
+  MakePair("OPEN", INC_OUTLET_TYPE::OPEN)
 };
 
 /*!
@@ -2493,6 +2532,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 . */
 };
 
 /*!
@@ -2524,6 +2564,9 @@ enum 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 communication. */
+  MASS_FLUX	       ,  /*!< \brief Primitive variable communication. */
+  MOM_COEFF,
 };
 
 /*!