# When to stop iterations in SOR solver for 3D Poisson equation

I'm writing a solver (in C) for 3D incompressible fluids, using the finite-differences method, and I'm finding a somewhat surprising behaviour: the solver provides "good-looking" solutions, but however pressure values keep increasing if you perform more iterations in the SOR solver.

Well, maybe this behaviour makes sense after all, because I'm applying Neumann BCs for pressure in all the boundaries (I'm doing it because that's the advice from the CFD books I'm studying), and so I guess that if you apply Neumann to pressure in all the boundaries, you end up having infinite solutions, so maybe always-increasing pressures is not as unexpected as I'm thinking (note: for velocities I'm applying Dirichlet in some boundaries and Neumann in others).

Because of this behaviour, I don't know when to stop the SOR iterations. In each iteration, I obtain the L2-norm of the pressure values after the iteration (let's call it p2norm) and the L2-norm of the residual after the iteration (let's call it res2norm). Then I obtain a dimensionless residual, dres=res2norm/p2norm.

For the model I'm testing (which is 25x25x100 cells), I found that it takes about 60 iterations for dres becoming smaller than 0.1, about 600 iterations for it becoming smaller than 0.01, about 6000 iterations for 0.001, and 60000 iterations for 0.0001, ...and certainly this looks a bit unexpected to my eyes.

The system never gets unstable (and the residual always decreases in every iteration, with no exception), but pressure values keep increasing. There's a difference, though, and it's that the more iterations I do, the more uniform the pressures get. For example, after 60 iterations, I get noticeable gradients in the pressure (min_pressure=4 to max_pressure=100). After 600 iterations, min_pressure=150 to max_pressure=400. After 6000 iterations, min_pressure=2300 to max_pressure=2800. After 60000 iterations, min_pressure=25300 to max_pressure=25700. As you can see, pressure tends to get uniform across the domain if I do a lot of iterations, although the range from the minimum to the maximum seems to stabilize around 400 or so.

I have tried with different values for the omega overrelaxation parameter (1, 1.2, 1.5, and 1.7), and the behaviour is always analogous (the orientative values I provided are from using 1.7, but the overall behaviour is similar in all the omegas I tried).

Do you have any suggestions on what other things I could check, or about something that I might be doing wrong? Have you ever seen a behaviour like this? In case it is correct, when should I stop the SOR iterations?