[WIP] Pressure Based Solver by thijsaalbers · Pull Request #2812 · su2code/SU2

thijsaalbers · 2026-05-07T13:05:07Z

Proposed Changes

The current work (part of GSoC) provides a working version of a pressure-based algorithm for the incompressible flow solver as an alternative to the existing Density-based solver. Below, the reader may find the algorithm which has been implemented, as well as the current progress of the code and challenges. All the way at the bottom one can find performance comparisons between the DB and PB solvers for some test cases.

Algorithm

A lot of versions of pressure-based algorithms exist, and many different versions can be implemented. Here, we opt for versions of the original SIMPLE/PISO algorithm, it is briefly defined here for clarity.

First, the momentum equations are solved, starting from the previous time step's velocity$\vec{u}^{(0)}$, pressure $p^{(0)}$, and face velocity $\vec{u}_f^{(0)}$. The resulting momentum is the predicted momentum, here its discretized form is shown, as its coefficients are used in the subsequent equations

$A_p\vec{u}_p^{(1)}+\sum_n A_n \vec{u}^{(1)}_n=-\frac{V}{\rho}\nabla p^{(0)}+S_m$

The subsequent momentum is not necessarily incompressible, the pressure correction equation can be derived by rewriting it as follows, using a term often called H by A

$\vec{u}_p^{(1)}=-\frac{\sum_n A_n \vec{u}_n^{(1)}}{A_p}-\frac{V}{\rho A_p}\nabla p+S_m=\frac{H(\vec{u}^{(1)})}{A}-\frac{V}{\rho A_p}\nabla p^{(0)}+S_m$

$\vec{u}_p^{(2)}=\frac{H(\vec{u}^{(1)})}{A}-\frac{V}{\rho A_p}\nabla p^{(1)}+S_m$

Note how a simplification is used here where the HbyA term is neglected. Subtracting these two equations yields the first pressure correction equation for $p'$ as

$\nabla \cdot \left( \frac{V}{\rho A_p}\nabla p'\right)=\nabla \cdot \vec{u}^{(1)}$

Make note that for the divergence here, we require the face mass fluxes, which are computed using Rhie-Chow interpolation to avoid odd-even decoupling. After the equation is solved, using the pressure correction $p'$, the pressure and momentum are corrected according to

$p^{(1)} = p^{(0)} + p' ,\quad u^{(2)} = u^{(1)}+u'.\quad u'= -\frac{V}{\rho A_p}\nabla p'$

So far, this is equal to a pseudo-transient version of the SIMPLE algorithm. This algorithm however suffers from a very tight stability condition on the time-step size. Therefore, multiple pressure corrections can be applied, which for two corrections is originally called the PISO algorithm.

The second pressure correction does not neglect the HbyA term, which then results in the equation

$\nabla \cdot \left( \frac{V}{\rho A_p}\nabla p'\right)=\nabla \cdot \vec{u}^{(2)}+\nabla\cdot\left(\frac{H(\vec{u}')}{A}\right)$

And the new correction equations are defined as

$p^{(2)} = p^{(1)} + p' ,\quad u^{(3)} = u^{(2)}+u'.\quad u'= \frac{H(\vec{u}')}{A}-\frac{V}{\rho A_p}\nabla p'$

Note that HbyA here uses the previous velocity correction and is thus the same quantity as the one used in the second pressure correction equation. Later pressure correction equations follow analogously.

Progress:

Pressure-based solver added as alternative to density-based solver for the incompressible flow equations.
The pressure-based solver has only been tested for constant density cases for basic Navier-Stokes and Euler flow.
The pressure-based solver is implemented based on a pseudo time-stepping approach to remain consistent with the other solvers in SU2. The pressure-based solver is currently set to the SIMPLE algorithm by default, with options for SIMPLEC and PISO available.
The Poisson solver is a major bottleneck in the computation speed. To account for this, an option was added to use a different linear solver and preconditioner for the poisson solver.
For details on the implementation of the algorithm and the responsibility distribution please see the file CPBFluidIteration.cpp.

Issues

Performance:

The Poisson solver can sometimes struggle a lot due to high Reynolds numbers and fine meshes, and thus require a ridiculous number of iterations to converge reasonably. Possible fixes include adding multigrid support or a DIC preconditioner (far less efficient). Multigrid support is tricky as SU2 currently only considers multigrid for the main (flow) solver and not for auxiliary solvers.
Convergence issues with RANS (SA and SST) on fine meshes with high Reynolds numbers. Tests have shown that cases such as flow over a flat plate converges fine. However, external aerodynamic cases such as the naca0012 RANS test case do not converge well at all. The convergence does slightly improve when we switch out the mesh for a more uniform unstructured mesh without large aspect ratio cells in the wake of the airfoil, although this only slightly helps. The flat plate turbulence test case also uses large aspect ratio cells so this is not the sole issue. The RANS solver also often requires many iterations of the Poisson solver to converge reasonably, this is however not the reason for the lack of convergence.
Periodic boundary conditions have not been implemented/tested at all as of yet.
Any code related to adjoints has not been considered at all either.

Code:

Parallelization with OMP gives wrong results, MPI however does work as expected.

TODO list

Fix issues mentioned above
Docs page

Related Work

This work is based on earlier attempts by Nitish Anand (2024) and Akshay Koodly (2021), see feature branches feature_PBFlow_V8 and feature_Pressure_based respectively. Also see PR #2210

PR Checklist

I am submitting my contribution to the develop branch.
My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
I have added a test case that demonstrates my contribution, if necessary.
I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

Result showcase

Inviscid Hydrofoil

Convergence history of the inviscid flow around a hydrofoil at a 5 degree aoa.

The pressure coefficient along the surface of the hydrofoil at a 5 degree aoa and the corresponding lift coefficients, X-FOIL predicts C_L=0.6.

Lid Driven Cavity

Convergence history of the lid driven cavity problem, note that CFL=60 is the highest stable CFL for the PB solver, whereas the DB solver does not have this CFL related stability issue.

Flatplate RANS

The skin friction coefficient for turbulent flow over a (rough) flat plate with SA.

github-advanced-security

CodeQL found more than 20 potential problems in the proposed changes. Check the Files changed tab for more details.

Implements the full class structure required for the pressure-based solver in a minimal form. The code compiles and runs but does not yet contain any numerical/physical implementation. Future work will focus on implementing solver logic. Note: In the previous attempts (see related work) there is noticeable code duplication between CIncEuler and CPBIncEuler. The final architecture may be revised depending on how the implementation of CPBIncEuler evolves.

Centered residual is not yet implemented as the old code did not have a working version. Upwind residual is functional but requires cleanup and move to more appropiate file.

Variables are a work in progress and still include some commented out code related to energy and pressure

The current viscous residual of the poisson equation uses an isotropic diffusion coefficient, this should not be the case and should thus be replaced

thijsaalbers · 2026-07-03T12:44:22Z

I will have to review the code myself first, it is not ready for review as of yet.

thijsaalbers

The current version of the solver is ready for review. It is however not yet finished and there is a list of known issues that I am still working through. I have updated the PR description with the current status and known limitations; please read that before reviewing the implementation.

pcarruscag

There are some errors in the implementation, which I made myself when I first implemented SIMPLE, and they are easy to make because certain simplifications look equivalent but are not.
For professional reasons I should discuss this only in terms of my prior work.
Here is an implementation that I think is mostly correct: https://github.com/pcarruscag/Flow-Solver-Experiments/blob/main/lib/flow/src/flow_simple.cpp
It's written in a funny way because it was my first exercise into porting something to GPUs, but follow the operations, and their order, not the code itself.
The main error here is that you do not carry nor correct the face velocities (or face mass flows), nor use them to discretize momentum and scalars.
You correct the nodal velocities (equivalent to cell centers) but that is not the same thing.
I was not nice enough to my future self to document why in my old code, but I'm somewhat confident an AI can explain the nuance.

thijsaalbers · 2026-08-03T17:58:19Z

@pcarruscag Thanks for the very usefull explanation of the issue. I think I now more or less fixed it. I also rewrote CPBFluidIteration.cpp, and I think the flow of the algorithm (and the changes I now made to it) should be clear in there.

…is still wrong

…ings

pcarruscag · 2026-08-21T03:54:32Z

+   * \param[in] poisson - If preconditioner should use settings defined for the poisson solver.
   */
  void Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn,
                  bool EdgeConnect, CGeometry* geometry, const CConfig* config, bool needTranspPtr = false,
-                  bool grad_mode = false, bool allow_quant = false);
+                  bool grad_mode = false, bool allow_quant = false, bool poisson = false);


Instead of grad_mode and poisson, we should just pass the enum value that we want to override.
Using a std::optional.

pcarruscag · 2026-08-21T03:54:47Z

  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.*/


Suggested change

Pressure_based; /*!< \brief FLag to check if we are using a pressure-based system.*/

Pressure_Based; /*!< \brief Flag to check if we are using a pressure-based system.*/

pcarruscag · 2026-08-21T03:56:32Z

+  /*!
+   * \brief Get the relaxation coefficient of the pressure correction for SIMPLE solver.
+   * \return relaxation coefficient of the pressure correction for SIMPLE solver
+   */
+  su2double GetRelaxation_Factor_Pressure(void) const { return Relaxation_Factor_Pressure; }
+
+  /*!
+   * \brief Get the relaxation coefficient of the momentum correction for SIMPLE solver.
+   * \return relaxation coefficient of the momentum correction for SIMPLE solver
+   */
+  su2double GetRelaxation_Factor_Momentum(void) const { return Relaxation_Factor_Momentum; }
+
+  /*!
+   * \brief Get the coefficient for the removal of the transient term in the poisson solver coefficients.
+   * \return coefficient for the removal of the transient term in the poisson solver coefficients.
+   */
+  su2double GetTransient_Term_Removal_Factor(void) const  { return Transient_Term_Removal_Factor; }
+
+  /*!
+   * \brief Verify if there is mixing plane interface specified from config file.
+   * \return boolean.
+   */
+  bool GetBoolAutomaticRelaxationFactors(void) const { return AutomaticRelaxationFactors; }


Bundle the SIMPLE-specific options in a struct like we do for turbulence models and ILU options please, that way there is a singe funtion to return them.

pcarruscag · 2026-08-21T04:01:23Z

  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. */
+  MOM_COEFF            ,  /*!< \brief Momentum coefficient for the pressure-based poisson solver. */


Suggested change

MOM_COEFF , /*!< \brief Momentum coefficient for the pressure-based poisson solver. */

MOM_COEFF , /*!< \brief Momentum coefficient for the Rhie-Chow scheme. */

pcarruscag · 2026-08-21T04:21:45Z

    /*--- Dirichlet condition for temperature at far-field (if energy is active). ---*/
+    if (!pressure_based) {

-    V_infty[prim_idx.Temperature()] = GetTemperature_Inf();
+      V_infty[prim_idx.Temperature()] = GetTemperature_Inf();


Same comment about just one if with continue in this case. Makes it a lot easier to review what is going on.

pcarruscag · 2026-08-21T04:22:57Z

+  /*--- Mass flux is computed over all edges ---*/
+  for (auto color : EdgeColoring) {
+    SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize))
+    for (auto k = 0ul; k < color.size; ++k) {
+      auto iEdge = color.indices[k];


You also don't need the coloring here since you write to the edges directly.
Coloring is for writing to the points of the edge.

pcarruscag · 2026-08-21T04:23:33Z

+      }
+      /*--- 1. Interpolate the pressure gradient based on node values ---*/
+      for (iDim = 0; iDim < nDim; iDim++) {
+        Grad_Avg = 0.5*(nodes->GetGradient_Primitive(iPoint,0,iDim) + nodes->GetGradient_Primitive(jPoint,0,iDim));


use the primitive indices class of the solver instead of hardcoding the index (0) please

pcarruscag · 2026-08-21T04:25:32Z

+      /*--- Rhie Chow interpolation ---*/
+      Coord_i = geometry->nodes->GetCoord(iPoint);
+      Coord_j = geometry->nodes->GetCoord(jPoint);
+      dist_ij_2 = 0.0;
+      for (iDim = 0; iDim < nDim; iDim++) {
+        Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim];
+        dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim];
+      }
+      /*--- 1. Interpolate the pressure gradient based on node values ---*/
+      for (iDim = 0; iDim < nDim; iDim++) {
+        Grad_Avg = 0.5*(nodes->GetGradient_Primitive(iPoint,0,iDim) + nodes->GetGradient_Primitive(jPoint,0,iDim));
+        GradP_in[iDim] = Grad_Avg;
+      }
+
+      /*--- 2. Compute pressure gradient at the face ---*/
+      /*--- Eq 15.62 F Moukalled, L Mangani M. Darwish OpenFOAM and uFVM book. ---*/
+      GradP_proj = 0.0;
+      for (iDim = 0; iDim < nDim; iDim++) {
+        GradP_proj += GradP_in[iDim]*Edge_Vector[iDim];
+      }
+      if (dist_ij_2 != 0.0) {
+        for (iDim = 0; iDim < nDim; iDim++) {
+          GradP_f[iDim] = GradP_in[iDim] - (GradP_proj - (nodes->GetPressure(jPoint) - nodes->GetPressure(iPoint)))*Edge_Vector[iDim]/ dist_ij_2;
+        }
+      }


Can you extract a common helper from the numerics classes? this is similar to the corrected gradient calculation.

pcarruscag · 2026-08-21T04:28:35Z

+
+    /*--- Velocity corrections ---*/
+    for (iDim = 0; iDim < nDim; ++iDim) {
+      nodes->SetSolution(iPoint, iDim + 1, nodes->GetSolution(iPoint,iDim + 1) + alpha_u * velocityCorrection[iPoint][iDim]);


This isn't how alpha_u is typically used in SIMPLE

Thats correct. Only as we currently use only pseudo-time stepping, there is not really a direct method for the user to influence the velocity corrections. It's not entirely necessary either indeed as the user can alter the CFL, which indirectly influences the time step size, which in turn indirectly functions as a form of alpha_u. So it is removable if you prefer, it doesn't matter that much.

We could also add an option to change the time step size based on the value of alpha_u and the jacobian matrix, which is mathematically equivalent to using a steady-state SIMPLE algorithm (with alpha_u as the underrelaxation value). This would however completely side-step the CFL parameter. So would this even be desired?

If CFL works ok and alpha doesn't do much, I would say remove alpha.
Traditional implementations do not relax the velocity and mass flow corrections directly.

Initial configuration options and a test case

64d4213

thijsaalbers self-assigned this May 7, 2026

thijsaalbers added PR: draft changelog:feature GSoC Google Summer of Code labels May 7, 2026

Numerics and driver additions, WIP

1a7e82a

github-advanced-security Bot found potential problems May 9, 2026

View reviewed changes

pcarruscag reviewed May 9, 2026

View reviewed changes

Comment thread SU2_CFD/src/numerics/poisson.cpp Outdated

Comment thread SU2_CFD/src/numerics/pbflow.cpp Outdated

Comment thread SU2_CFD/src/numerics/pbflow.cpp Outdated

thijsaalbers added 2 commits May 12, 2026 13:24

Remove outdated numerics

8d0bcf9

github-advanced-security Bot found potential problems May 12, 2026

View reviewed changes

thijsaalbers added 2 commits June 3, 2026 11:30

Add numerics

a1805d4

Centered residual is not yet implemented as the old code did not have a working version. Upwind residual is functional but requires cleanup and move to more appropiate file.

Add variables and output

d28be3b

Variables are a work in progress and still include some commented out code related to energy and pressure

thijsaalbers closed this Jun 4, 2026

thijsaalbers added 12 commits June 10, 2026 14:16

MPI quantities

0b7f37b

Added momentum equations

e542979

First working version of pressure correction

ef875fb

The current viscous residual of the poisson equation uses an isotropic diffusion coefficient, this should not be the case and should thus be replaced

Added PISO and linear solver options

1931fa9

Altered some config parameters

81f47a3

Removed duplicate code in iteration

30933b8

Small fix for poisson preconditioners

425c39d

Small cleanup on numerics

51827d8

forgotten meson file

c43dbba

Workaround for symmetry BC

c7aadfa

Latest working version of the solver

e2ec0e2

small change

d588424

thijsaalbers reopened this Jul 3, 2026

thijsaalbers added 2 commits July 3, 2026 14:51

Merge branch 'develop' into feature_pressure_based_2026

9207743

Fixed loading of restart file

9ca3e00

thijsaalbers added 2 commits July 20, 2026 23:52

fix compile issue with OMP

a8d8a8d

fixed compile issue with OMP

5ed114c

thijsaalbers commented Jul 21, 2026

View reviewed changes

thijsaalbers requested a review from pcarruscag July 21, 2026 12:29

Small cleanup and testcases

f3536fd

thijsaalbers marked this pull request as ready for review July 24, 2026 14:21

pcarruscag reviewed Jul 29, 2026

View reviewed changes

thijsaalbers added 4 commits July 31, 2026 10:51

fixed tests

a992388

Fix for euler solver issue

3c88ce0

Refactoring to make use of edge velocities

4741608

bounded scalar, and small fixes

909870e

thijsaalbers added 17 commits August 8, 2026 14:58

fixed PISO algorithm which fixes the stability issues at higher CFL

2999a8b

Merge branch 'develop' into feature_pressure_based_2026

be74da7

Fixed conflicts

79fbf38

Added unsteady density to pressure based solver and config template

3ad2d75

Unified DB and PB solvers, density variable for PB remains a WIP

f85b849

Fixed order of variables for thermo purposes

f67a94c

minor style changes

b5ba239

Fixed MPI for PISO algorithm

17aef33

Improved issues with OMP, still doesnt work completely

2ab2c9b

Fixed non-unity density

3810647

Updated test cases

80fdc42

Merge branch 'develop' into feature_pressure_based_2026

1acdc05

Refactored solution vector of the solver to match DB solver, density …

c2e50fd

…is still wrong

Fixed density for new solution vector and small clean up for surround…

0a8c7d2

…ings

fixed the OMP issue

85a310c

Added more, and more clear, configuration options

9bdbe57

Merge branch 'develop' into feature_pressure_based_2026

c29ad00

pcarruscag reviewed Aug 21, 2026

View reviewed changes

	Pressure_based; /!< \brief FLag to check if we are using a pressure-based system./
	Pressure_Based; /!< \brief Flag to check if we are using a pressure-based system./

	MOM_COEFF , /!< \brief Momentum coefficient for the pressure-based poisson solver. /
	MOM_COEFF , /!< \brief Momentum coefficient for the Rhie-Chow scheme. /

Uh oh!

Conversation

thijsaalbers commented May 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed Changes

Algorithm

Progress:

Issues

Performance:

Code:

TODO list

Related Work

PR Checklist

Result showcase

Inviscid Hydrofoil

Lid Driven Cavity

Flatplate RANS

Uh oh!

github-advanced-security AI left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

thijsaalbers commented Jul 3, 2026

Uh oh!

thijsaalbers left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

thijsaalbers commented Aug 3, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

thijsaalbers commented May 7, 2026 •

edited

Loading