This post describes a brief overview of a small investigation I did of realtime fluid simulation for interactive purposes, written in C++. The project implements the Smoothed Particle Hydrodynamics model of fluid dynamics described by Matthias Muller et al (2003). Fluids are simulated and rendered in real-time, making them appropriate for use in interactive applications such as video games or physics sandboxes.

Demo

A short video demonstrating the fluid, boundary conditions, and rendering using marching squares or the debug render (density):

Physics introduction

Smoothed Particle Hydrodynamics is a Lagrangian method for simulating the movement and interactions of fluids. Unlike Eulerian grid-based approaches, SPH discretises a fluid into distinct particles.

There are five three distinct phases for each step of the simulation:

1. The density of each particle is evaluated by summing the contribution of each nearby particle - the more particles nearby, the greater the density of that particle (note: this simulation is therefore clearly compressible). The contribution of each particle is smoothed by a smoothing function - which weights the contribution of other particles based on their distance from the current particle.
2. The pressure force is calculated for each particle by again summing the contribution of each nearby particle - this time, as a vector. Particles will push away other particles. Denser particles are less able to push away their neighbours.
3. The viscosity force is calculated for each particle, based on the velocity difference between the particle and it's neighbours.
4. The particles are integrated through time using Euler's method (less accurate than other methods, but the simplest to calculate).
5. Boundary collisions are resolved by pushing the particles out of the wall.

Implementation Features

Naively, to calculate the density of a particular particle would require summing against all other particles, an O(n^2) operation on all particles. This time complexity is prohibitive for real time applications using thousands of particles.

The smoothing function is only non zero when the particle is within the smoothing radius, h, of other particles. This lends itself well to optimisation. I have integrated a grid-based to ensure that when computing the density, pressure and net-force on a particle, only the particles immediate neighbours are considered. This grid is designed to trade off storage in favour of cache-speed, storing copies of particles in a vector rather than references to these particles. In conjuction with the effective limit of how many particles can physically inhabit the grid, the grid linearises the time complexity of the calculations.

Particles in SPH are discrete. In order to display a contiguous surface, a mesh must be generated. I have implemented a marching squares algorithm, which produces a mesh by evaluating the density field of the fluid at discrete points. The actual generation of the mesh was performed in an OpenGL geometry shader, which uses the parallel computing power of the GPU to add vertices to a surface. This allows for generating high resolution surfaces.

The use of a grid also makes for efficient multithreading. The grid cells can be divided evenly by worker threads, which each are responsible for calculating certain particles. This is used to great effect in the renderer, which uses multiple threads to evaluate the fluid density.

Modern C++

I've made a conscious effort throughout this project to use modern C++ features. These include:

• RAII
• Smart pointers to avoid memory leaks
• Smart iterators to avoid bounds errors
• Platform independent threads and mutexes

Future Work

The main future work to be explored in this project is to make use of the GPU, either in CUDA, Vulkan, or a Compute Shader, to perform all fluid calculations on the GPU. That would greatly accelerate the computation, and allow orders of magnitude more particles. In addition, although in theory the simulation code would work in three dimensions, the rendered requires extension to support 3D surfaces.