[andreasp]

Hi everyone!

For demonstration purposes, I have implemented a little 2D incompressible flow solver using the SIMPLE(R) algorithm by Patankar on a staggered grid. My serial solver version seems to work fine, I can simulate lid-driven cavity and different channel flows.

Now my next step is to parallelize the code using a fairly simple OpenMP approach. So the major parts to parallelize are a couple of loops over all the grid cells and faces (purely cartesian mesh). For example, in the serial code version I have one loop over all cells to compute the pressure correction. The computation of pressure correction for each cell involves the pressure correction values of its four neighbor cells.

for all cells i:
    press_corr[i] = func( press_corr[west], press_corr[east], press_corr[south], press_corr[north] )

Typically, two of these neighbor values are the "old" values from last iteration, and the other two are already the new values (in those cells which have already been visited in the same loop). Parallelizing such a loop without care is not necessarily dangerous (because each value is modified by only one thread), but nevertheless can lead to non-deterministic results, because it is not clear in which order the cells are visited.

To circumvent this, I decided to introduce an intermediate storage which keeps the new values. These new values are now computed exclusively from old values, and then the new values are copied back to the actual pressure correction array in a second loop:

for all cells i:
    press_corr_intermediate[i] = func( press_corr[west], press_corr[east], press_corr[south], press_corr[north] )

for all cells i:
    press_corr[i] = press_corr_intermediate[i]

These two for-loops can be parallelized trivially without generating any non-deterministic results. So far, so good. However, here comes my problem:
With this modification, the numerical method does not converge anymore! My residuals blow up.
As long as the pressure correction calculation depends partially on new (current iteration) neighbor values, the method converges. But as soon as I use exclusively old neighbor values (from last iteration), it does not converge anymore. (This is completely independent of any parallelization! Just the same happens in serial execution.)

Can anyone give me a hint on what might be going on here?

Many thanks in advance!

Cheers,
Andreas

[mprinkey]

The serial version is doing Gauss Seidel iteration. The two-loop parallel version is doing Jacobi iteration. These are not the same iterative scheme.

If you don't have a fixed pressure boundary condition, the pressure (correction) system will be formally indeterminate as I've discussed here a few months ago. GS iteration may converge while Jacobi may not. It might require you to fix one pressure value to make the system determined.

A better option is to do red/black Gauss Seidel. That avoids the parallel dependency issue and should give convergence more similar to the original.

[andreasp]

Thank you very much for your quick answer!
That totally makes sense. Indeed, the system that is not solvable is a lid-driven cavity problem which does not have any fixed pressure values.

I guess I will take a look at the red/black GS scheme.

Thanks again!

[Ramkumar21194]

Hi, i got similar problm while solving 2d viscid flow using simple algorithm, but the solution diverges and the values of pressure and temperature reaches higher values till infinity... plz suggest me on this anyone...the higher value arises from pressure correction equations...
Thanks in advance.