Eulerian incompressible smoothed particle hydrodynamics on multiple GPUs

Recent advances in the development of Eulerian incompressible smoothed particle hydrodynamics (EISPH), such as high-order convergence and natural coupling with Lagrangian formulations, demonstrate its potential as a meshless alternative to traditional computational fluid dynamics (CFD) methods. This work aims to address one of the major outstanding limitations of EISPH, its relatively high computational cost, by providing an implementation that can be deployed on multiple graphics processing units (GPUs). To this end, a pre-existing multi-GPU version of the open-source Lagrangian weakly-compressible code DualSPHysics is converted to an EISPH formulation and integrated with an open-source multi-GPU multigrid solver (AmgX) to treat the pressure Poisson equation implicitly. The integration of AmgX within DualSPHysics presents a significant challenge, since AmgX is designed for distributed systems and therefore conflicts with the single-node shared memory design of the multi-GPU DualSPHysics code. The present implementation is validated against well-known test cases, showing excellent agreement with benchmark solutions and demonstrating second-order convergence. A detailed profiling and performance testing is also presented to investigate memory consumption and scaling characteristics. The results show approximately 87%–95% strong scaling efficiency and 92%–94% weak scaling efficiency in both 2D and 3D on up to four GPUs. Large spikes in memory consumption during the initialisation of the linear solver library are found to impede full utilisation of the device memory. Nevertheless, the present implementation is shown to permit problem sizes on the order of 69.3 million (2D) and 20.8 million (3D) particles on four GPUs (128 GB total device memory), which is beyond what has previously been reported for incompressible SPH on GPUs, where the PPE is treated implicitly, and demonstrates its potential as an alternative to traditional CFD methods.