[panagiat]

Hello!

I recently made a steady-state FVM Incompressible 2D Navier Stokes solver using the SIMPLE algorithm for a collocated cartesian orthogonal grid in C++ by following the lectures by Sandip Mazumder

For solving the momentum/pressure correction equations, I am using the Gauss Seidel algorithm.
As for the residual calculation, I am using the L^2 norm, adding the error of each node and at the end taking the square root

As far as my understanding of the SIMPLE algorithm goes this is the general structure:

while u, v, p are not converged (outer iteration):
    Solve X momentum (with inner iterations)
    Solve Y momentum (with inner iterations)
    Solve pressure correction equation
    Correct stuff

I am going to use the terms initial/final residuals which are explained here

I am calculating the initial residuals as follows:
- Velocity: After the 1st inner iteration of the Gauss seidel algorithm of the momentum equations.
- Pressure: Calculating the mass imbalance after calculating the pressure correction coefficients right before solving the equation.

I am calculating the final residuals as follows:
- Momentum equations: The residual values of the last iteration of the Gauss Seidel algorithm.
- Pressure correction equation: The residual values of the last iteration of the Gauss Seidel algorithm.

At the moment I am only doing 1 inner iteration for each equation and the code converges for the lid driven cavity problem I am solving. Though it seems that nobody is doing it this way seeing diagrams of residuals commercial products create.
(maybe because this way the coefficients are calculated more times than needed?)
Thus, I decided to start adding more inner iterations and this is where trouble began.

I am solving the momentum equations by reducing the initial residual I calculated on the first inner iteration by an amount of orders (set at 2 orders of magnitude).
As for the pressure correction equation, I am starting with 0 pressure correction everywhere and I am solving with a tolerance of 4 orders of magnitude.

The momentum equations seem to be converging, but something weird happens with the pressure correction equation. Instead of the pressure correction values converging, it seems that the difference between each inner iteration is converging while the correction values keep on getting larger. This is obviously not the expected behaviour. On a closer look, the first outer iteration goes by because it just happens that the difference is within the tolerance, but on the second outer iteration it is pretty clear that something's wrong (which is probably why I got away with just 1 inner iteration).

Have I understood something wrong or is it simply an error in the way I am solving the pressure equation (or the coefficients themselves)?
Any ideas?