Try: Esc Space r Left Mouse Right Mouse

Three weeks ago I saw a video where Sebastian Lague, a Unity developer and YouTuber, created an interactive fluid simulation in Unity. The project hooked me and I decided to make my own in C++. Here’s how it went.

How to simulate a fluid

In fluid simulation, there are two big classes of simulation:

  1. Eulerian
  2. Lagrangian

Either of these models can be used to simulate either liquid or gas. My simulation attempts to simulate liquid meaning that it should have a more or less constant density to achieve incompressible behavior. I saw the term “semi-incompressible” used in a few papers which is how I’ll justify the clearly variant density of the fluid above.

Eulerian simulations are grid-based simulations that compute inflow and outflow per cell to get incompressible behavior. These simulations tend to be highly accurate and useful for scientific or engineering applications.

Lagrangian simulations are particle-based and rely on tracking individual pieces of fluid. These simulations tend to be used in rendering due to their ability to simulate fluid in dynamic environments.

The technique I used is a widely adopted lagrangian method called Smoothed Particle Hydrodynamics (SPH). This technique applies a “smoothing kernel” to estimate fluid properties at some location. Concretely, a smoothing kernel may estimate density by summing the masses of the particles that fall within the “smoothing radius” and weighting those masses by the scalar returned by the smoothing kernel.

This diagram illustrates the concept well.

Source: Taken from [1]

Source: Taken from [1]

Using SPH we can compute a density and a density gradient then offset the density gradient by a target density. Now all we have to do is move each particle in the direction of the density gradient and the fluid should converge to the target density. We use a constant pressure value to determine how strongly to push each particle in the direction of the density gradient.

In practice, this can be a very finicky process, where small changes to the number of particles, target density, or pressure multiplier can lead to an unstable simulation. I have yet to find a way to systematically “tune” fluid behavior. You can play with these values from the debug menu by pressing Esc in the example above and don’t forget that you can reset the values by pressing r.

Optimization

After implementing a basic SPH simulation I couldn’t help but notice that my simulation wasn’t very fast. Framerate took a serious nosedive if anymore than 400 particles were drawn. It’s not hard to guess that single-threaded, un-optimized simulation won’t be very fast but it’s good to do more than guess so I cracked open a profiler to see where my program was spending all its time.

MacOS Instruments Profile

MacOS Instruments Profile

Sure enough, over 90 percent of time was being spent in the simulation step function which is doing all the simulation math.

Fortunately, there are some low-hanging fruit for increasing the efficiency in this step. The first and most obvious optimization is to figure out how to only do math on particles that might influence each other.

I’m not the first person to have this problem and the common solution is to create a grid where each cell tracks the particles it contains. If we set the grid cell size to be the same size as the smoothing radius then we will never need to check more than nine cells for any particle.

Source: Grid comparison diagram from [2]

Source: Grid comparison diagram from [2]

If you imagine the red particle in a corner of the center cell it is easy to see that the smoothing radius (denoted with h) would never reach outside the highlighted 9-cell region.

Implementing this change produced a massive speedup and allowed me to jump from several hundred particles to just under ten thousand. Unfortunately, I didn’t take a profile at this point so you’ll just have to take my word for it.

You can see the particles that are in this 9-cell region by using the “Adjacency view” from the debug menu. This will show particles in cells adjacent to your mouse. This was implemented using a C++ forward_iterator which lets me iterator over adjacent particles in a range-based for loop like this:

1
2
3
4
5
  void Fluid::computeDensity(const Vec &p) { 
    for (int i : this->grid.adj(p)){ 
      // do stuff
    }
  }

Parallelization

Though the grid optimization allows for far more particles, it’s not nearly fast enough to support a 3D simulation of any real size. For this, we’d need another order of magnitude speedup. The obvious candidate is to use the graphics card to compute how each particle interacts with its neighbors in parallel. This would effectively take the computation that requires n_particles * n_interactions_per_particle iterations and transform it into n_interactions_per_particle iterations that all happen at once.

A great option for doing this kind of parallel computation is CUDA. CUDA is a proprietary technology from Nvidia that provides a nice API for doing parallel GPU computation. Unfortunately, CUDA is only available on Nvidia graphics cards and my poor Macbook Air didn’t come equipped with this kind of hardware. Even if it did CUDA isn’t supported by web assembly.

Web assembly also foiled my next thought which was to use OpenGL 4.3’s compute shader feature. This feature gives somewhat analogous features to CUDA but through an open standard that is more widely supported.

This leaves me with using the old OpenGL vertex and fragment shaders to do particle math. My brain doesn’t yet support doing this implementation. Though I did find an excellent post by nicebyte demonstrating how a similar project can be done and maybe at some point in the future, I’ll learn OpenGL.

C++ and Web Assembly

To simulate a fluid we need to render it somewhere. There are lots of options for drawing in C++. I reached for SFML for its ease of use and my lack of experience with OpenGL.

SFML provides an easy API with nice function like drawCircle which I used to render each particle. This worked great for getting the math working. It was also easy to setup Imgui as a debug menu.

When I decided to use Emscripten to compile to web assembly I discovered that SDL (a different render context manager) is supported by while SFML is not. The switch to SDL was relatively painless, though I did forgo the standard OpenGL graphics tooling to simplify the transition.

Emscripten is a full compiler toolchain that has been used to port lots of old C/C++ projects to web assembly. In fact, it can be used to webify most projects that can be compiled by Clang / LLVM.

As C++ tools go Emscripten is easy to use, though it does have a mountain of compiler flags to read through and I’m still not entirely sure I’ve picked the best (or even correct) settings for my project.

Check out the source code on my Github or a full-screen version of the demo here.

Citations

[1] “Smoothed-particle hydrodynamics,” Wikipedia, https://en.wikipedia.org/wiki/Smoothed-particle_hydrodynamics

[2] D. Koschier, J. Bender, B. Solenthaler, and M. Teschner, “Smoothed particle hydrodynamics techniques for the physics based simulation of fluids and solids,” Eurographics DL Home, https://diglib.eg.org/items/e6d5e2df-9b18-4a1c-9f35-35c14412200b