Stable Fluids on CPU (macOS, OpenMP, OpenGL 4.1)

I’ve always been interested in the mathematics behind shaders and Computer Graphics, so when I stumbled upon Pavel Dobryakov’s WebGL fluid simulation (https://paveldogreat.github.io/WebGL-Fluid-Simulation/) I knew I had to try to create something like it.

I tried to build this fluid simulator for regular, consumer-grade CPUs (such as the one on my MacBook Air), however, since OpenGL on MacOS stops at 4.1 (no compute shaders), and I didn’t feel like learning Metal for this, I decided to use OpenMP to see how close I could get using only the CPU.

---

Dependencies

The build script fetches most UI and windowing dependencies automatically during the build process.

Component	Purpose
CMake (3.20+)	Orchestrates the build and manages FetchContent modules.
OpenGL 4.1	Final supported legacy OpenGL version on macOS; used for hardware-accelerated texture blitting and ImGui rendering.
Git	Required by CMake to clone dependency repositories.
OpenMP (optional)	CPU parallelization for the simulation. Improves performance at higher grid sizes, but the code runs well without it at lower settings. Install via Homebrew: brew install libomp.

Bundled via CMake FetchContent

Library	Role
GLFW	Window creation, context management, and input handling.
ImGui	Real-time UI for parameter tuning and simulation control.

Build & Run (macOS)

The code includes my build.sh script that automates the CMake configuration. It's specifically tuned for macOS, and handling the search for libomp (if installed via brew) and setting the appropriate Clang compiler flags for OpenMP support.

git clone https://github.com/AndrewBlackwell/fluid-sim.git
cd fluid-sim
chmod +x build.sh
./build.sh
./build/fluid

Note: The script logs most output to build/logs/ to keep your terminal clean.

---

Underlying Maths (for nerds)

This solver implements the incompressible Navier-Stokes equations for homogeneous fluids. The solver follows the stable, grid-based approach described by Jos Stam (“Stable Fluids”, SIGGRAPH 1999).

The evolution of the velocity field over time is described by the following, subject to the incompressibility constraint:

$$\frac{\partial \mathbf{u}}{\partial t} = \underbrace{-(\mathbf{u} \cdot \nabla)\mathbf{u}}_{\text{Advection}} + \underbrace{\nu \nabla^2 \mathbf{u}}_{\text{Diffusion}} - \underbrace{\nabla p}_{\text{Pressure}} + \underbrace{\mathbf{f}}_{\text{Ext. Forces}}$$

In src/fluid.cpp, the Fluid2D::step function breaks this equation down into sequential operations using Operator Splitting.

Advection (Semi-Lagrangian)

The non-linear term $-(\mathbf{u} \cdot \nabla)\mathbf{u}$ represents self-advection (the fluid moving itself). Standard Eulerian integration here is unstable, so I use a Semi-Lagrangian scheme, loosely as follows:

1. For every cell center $\mathbf{x}$, backtrace the position of the particle to where it was $\Delta t$ ago (dt in the code): $\mathbf{x}_{prev} = \mathbf{x} - \mathbf{u}(\mathbf{x}) \Delta t$.
2. Then, bilinearly interpolate the old velocity field at $\mathbf{x}_{prev}$ and assign it to the new grid cell.
3. See Fluid2D::advect; this method is unconditionally stable, allowing for larger time steps (dt) without the simulation blowing up.

Diffusion (Implicit)

The term $\nu \nabla^2 \mathbf{u}$ describes how fluid viscosity dampens velocity differences. An explicit solution would require tiny time steps to remain stable.

My code aims to solve this implicitly:

1. ($\mathbf{I} - \nu \Delta t \nabla^2$) $\mathbf{u}_{new} = \mathbf{u}_{current}$, resulting in a sparse linear system (a Poisson equation). This is implemented in Fluid2D::diffuse.
2. Solve the system using Gauss-Seidel relaxation inside Fluid2D::lin_solve.

Projection (Mass Conservation)

After advection and diffusion, the velocity field often becomes divergent (compressible). To enforce $\nabla \cdot \mathbf{u} = 0$, gotta perform a Helmholtz-Hodge Decomposition. This is handled in Fluid2D::project, which forces the pretty swirling vortices typical of incompressible flowing fluids.