Parallel simulations using MPI.jl

[marinlauber]

This pull request is a work in progress.

I open it now to see how we can add parallel capabilities to WaterLily efficiently, keeping in mind that ultimately we want to be able to do multi-CPU/GPU simulations.

I have run most of the examples/TwoD files with the double ghost again, and it works in serial.

There are a lot of files changed/added in this pull request and I will briefly describe the changes done.

changed files

• src/WaterLily: enables passing the type of the Poisson solver to the simulations (mainly to simplify my testing), preliminary MPI extension (not used, to be discussed).
• src/Flow.jl: implement double ghost cells and remove the special QUICK/CD scheme on the boundaries, as these are no longer needed.
• src/MultiLevelPoisson.jl: implement downsampling for double ghost arrays and change all utilities functions accordingly. Explicitly define the dot product functions to overload these with the MPI functions later on. Also, change the solver! function as the PoissonMPI.jl test did not converge properly with these Linfty criteria.
• src/Poisson.jl: add a perBC! call in Jacobi (not needed, I think) and adjust the solver.
• src/util.jl: adjust all the inside functions and loc to account for the double ghost cells. Adjust BC!, perBC! and exitBC! functions for the double ghost cells.

New files

• WaterLilyMPI.jl: contains all the function overload needed to perform parallel WaterLily simulations. Define an MPIGrid type that stores information about the decomposition (global for now) and the mpi_swap function that performs the message passing together with some MPI utils.
• MPIArray.jl: a custom Array type that also allocates send and receive buffers to avoid allocating them at every mpi_swap call. This is an idea for the final implementation and has not been tested yet.
• FlowSolverMPI.jl: tests for some critical part of the flow solver, from sdf measures to sim_step. Use with vis_mpiwaterlily.jl to see plot the results on the different ranks.
• PoisonMPI.jl: parallel Poisson solver test on an analytical solution. Use with vis_mpiwaterlily.jl to see plot the results on the different ranks.
• diffusion_2D_mpi.jl: initial test of MPI function, deprecated
• vis_diffusion.jl: use to the the results of diffusion_2D_mpi.jl deprecated
• ext/WaterLilyMPIExt.jl: my initial try of adding MPI as an extension, not used currently.
• test/poisson.jl a simple Poisson test, will be removed
• test/test_mpi.jl initial MPI test should be changed to use this instead.

The things that remain to do

• decide how to incorporate this into the main solver; this introduces many changes, like double ghosts, etc. We might want to consider whether this is another version or another branch or if we merge it into the main solver.
• try to minimize global operations in the code (for example, the AllReduce in Poisson.residuals!)
• proper debug (sometimes it fails with an MPI error; I haven't managed to track it yet, but I suspect it was from using @views() for the send-receive buffer. This could be avoided if we allocate the send and receive buffer with the arrays using something similar to what it is in the file MPIArray.jl
• Update the VTK extension to enable the writing of parallel files.
• Think about a GPU version of this...
• See how we can hide the MPI communication with array core computation to achieve perfect weak scaling (a la ImplicitGlobalGrid.jl)
• benchmark it and write proper tests

Some of the results from FlowSolverMPI.jl

basic rank and sdf check

zeroth kernel moment vector with and without halos

full sim_step check

[marinlauber]

@weymouth @b-fg @TzuYaoHuang My initial idea is to incorporate this into the main solver as an extension and use the custom type MPIArray and type dispatch to use the parallel version of the functions that need changing (all the functions in WaterLilyMPI.jl essentially). For example, ext/WaterLilyMPIExt.jl would define

BC!(a::MPIArray,A,saveexit=false,predir=())
    ...
end

The flow will be constructed using Simulation(...,mem=MPIArray), which will trigger the call of the correct functions. This will work for most of the functions, except loc, which might be annoying. Using the custom type will also mean that we won't allocate anything during MPI exchanges; the send and receive buffers are allocated as part of these arrays, so they can be accessed in the function that needs a transfer. We allocate only 2 send and receive buffers, with sufficient size to store any of the buffers in 3D.

We could also bind the MPIGrid to these array subtypes to avoid a global definition of the grid, but then I don't know how to deal with functions like loc.

[marinlauber]

this can be done upon initializing the writer to avoid global operation at every write.

[b-fg]

just FYI --project=. and just --project are equivalent.