Incompressible Fluid Simulation Parallelization with openMP and CUDA
Xuan Jiang, laurence lu, Linyue Song · 2022
We note that we base our initial serial implementation on the original code presented in Jos Stam’s paper and try several ways of parallelization of the computation. From this initial implementation, it was easiest to implement OpenMP. Because of the grid-based nature of the solver implementation and the shared-memory nature of OpenMP, the serial implementation did not require the management of mutexes or otherwise any data locks, and the pragmas could be inserted without inducing data races in the code. We also note that due to the Gauss-Seidel method, which in solving a linear system only requires intermediate steps, it is possible to introduce errors that cascade due to relying on neighboring cells which have already been updated. However, this issue is avoidable by looping over every cell in two passes such that each pass constitutes a disjoint checkerboard pattern. Because all dependencies are strictly in cardinal directions (as opposed to diagonal), this avoids relying on partially completed computations and guarantees correctness.While the serial implementation has highly parallelizable components, the components that are easily parallelizable are often interleaved with highly sequential components. To be specific, the set_bnd function for enforcing boundary conditions has two main parts, enforcing the edges and the corners respectively. The edges are easily parallelizable, but because there are only four corner cells, this is most easily done sequentially. However, this imposes a strange implementation where we dedicate exactly a single block and a single thread to an additional kernel that resolves the corners, but it’s almost not impacting the performance at all and the most time consuming parts of our implementation are cudaMalloc and cudaMemcpy. The only synchronization primitive that this code uses is __syncthreads(). We carefully avoided using atomic operations which will be pretty expensive, 1 but we need __syncthreads() during the end of diffuse, project and advect because we reset the boundaries of the fluid every time after diffusing and advecting. This is a GPU kernel level synchronization that should ensure this correctness. Otherwise, the code may decide to set boundary without finishing diffusing or advecting. We also had to change the serial code and include set_bnd_finish to finish the setting boundaries.