Extending the Simple Fluid Simulation

In the first post introducing The Construct, I showed how the library could be used to create a simple fluid simulation with little effort. The visual effect wasn’t really that awesome for a couple of reasons. We’ll extend that same “project” in this post and show how we can implement a couple of ideas and papers from Computer Graphics literature with little additional effort.

Artificial Viscosity

One of the main problems with the Semi-Lagrangian advection (SLA) approach is that, being only an approximate solution, it introduces error into the system. When we advect quantities that lie in a grid around using SLA, quantities are interpolated from nearby points and this leads to a kind of diffusion of the velocity, density, or whatever happens to be advected. This turns out to be a difficult problem to solve in general and there isn’t really any easy and efficient silver bullet to fix it.When the velocity field is advected using SLA, it causes a loss of energy and vorticity (which is a measure of spinning and rotation in the fluid). Since this spinning motion is the lively thing that we would like to capture in the computer graphics setting, we would like a way to mitigate this loss of vorticity.

One technique, called Vorticity Confinement, seeks to mitigate this effect by essentially measuring where vorticity is being lost, and to inject rotational motion back into the flow. While it’s not a perfect solution, it is pretty easy to implement and does give a more lively flow. Mathematically, the procedure of computing Vorticity Confinement is also pretty straightforward.

We first calculate the vorticity by taking the curl of the velocity

$\B \omega = \nabla \times \B u$

In The Construct, this is simply computed as

VectorField omega = curl(velocity);

Once we have the curl, we compute the the gradient of the magnitude of this vorticity, which will be a vector field pointing towards locations of high vorticity. Mathematically we have

$\B \eta = \nabla |\B \omega| $

$\B N = \B \eta / |\B \eta|$

which is handily translated to

omega = writeToGrid(omega, constant(Vec3(0,0,0)), domain);
VectorField eta = grad(length(omega));
VectorField N = eta / (constant(.00001) + length(eta));

Just a couple of notes here. First, we had to write omega to a grid before we created eta, for the simple reason that $\B \eta$ involves a gradient of the curl, which would require analytically taking two derivatives of the velocity, which isn’t possible in The Construct. By writing it to a grid, we can take another derivative without worry. Secondly, when we compute N, we had to add a small epsilon to the denomintator since it is possible (and even likely) that some places will have zero velocity/curl which would create a divide by zero. This just ensures that we do not divide by zero when we normalize this vector field.

Finally, the actual Vorticity Confinment force pushes perpendicular to this gradient like a paddle wheel churning up turbulence:

$\B f_{VC} = \epsilon h (\B N \times \B \omega)$

where $\epsilon$ is a user parameter adjusting the strength of the confinement forcing, and $h$ is the side length of a voxel in the grid. This means that the entire Vorticity Confinement force comes down to just a few more extra lines of code, which for the most part, mirror the mathematics almost exactly as described in the paper:

VectorField VorticityConfinement(
 VectorField velocity, Domain domain, float epsilon) { 
   VectorField omega = curl(velocity);
   omega = writeToGrid(omega, constant(Vec3(0,0,0)), domain);
   VectorField eta = grad(length(omega));
   VectorField N = eta / (constant(.00001f) + length(eta));
   return constant(epsilon * domain.H[0])) * cross(N,omega);

Applying this force to the original simulation gives us these results:


The result is pretty subtle, but if you looks closely, you will notice that there are finer features, and more importantly that the motion of the fluid does not dampen as quickly as the vanilla simulation.

Rayleigh-Taylor Instability

One interesting scenario that we can compute demonstrates a phenomenon called the Rayleigh-Taylor Instability. This instability refers to a shift to chaotic motion after a symmetry is broken across an interface of two varying densities. Essentially you can imagine a “heavy” material resting in a perfectly balanced way on top of a less dense material. So long as there is no net force across the boundary, the system will remain in a quasi-stable state. The moment a slight disturbance destroys this perfect symmetry, there is a net force which causes the motion to become chaotic.

We can mimic this scenario by create two densities, a red one on top that is heavy, and a blue one on the bottom which is lighter.

auto red = mask(dot(identity(), constant(Vec3(0,1,0))));
auto blue = mask(dot(identity(), constant(Vec3(0,-1,0))));

We only need to augment the force now so that red pushes downward, and blue pushes upward, and have each of these advect along the velocity field independently, since we have two different densities now.

for(...) {
  red = advect(red, velocity, delta_t);
  red = writeToGrid(red, constant(0.f), domain);
  blue = advect(blue, velocity, delta_t);
  blue = writeToGrid(blue, constant(0.f), domain);

  VectorField force = (blue - red) * constant(Vec3(0,1,0));

One small extra thing we can do is add a small bit of asymmetry to one of the densities so that the turbulent transition will occur more quickly. After red is initialized, we can just add a small bit of extra density to break things up:

VectorField x = identity(), N = constant(Vec3(0,0,1));
ScalarField extra = mask(constant(.05f) - length(x-N*dot(x,N)));
red = red + extra;

Putting these things together produces the following result:

Gridless Advection

This next example is rather curious in that we will only need to remove code in order to see its effect. Strictly speaking, we do not need to write the density into a grid at each time step. We have been writing the density into a grid at each time step because this means that the subsequent time step will only need to read from a grid to get a value for the density. If we do not write any of the density to a grid at each time step, we will construct a larger and larger expression. Writing density to a grid makes things fast, but we do pay a price. Because later reads from that grid require interpolation, diffusion takes place, which makes the density “smear” out as it flows. If we avoid writing density to a grid, then this diffusion is removed!

Removing those lines is easy, but be prepared for long render times. Two examples of the results of this are below, rendered from an angle with a version of my own renderer (approximately 2 hours to render the higher-resolution video):

SELMA and the Characteristic Map

Gridless Advection was really quite nice visually, but takes forever to compute, because each evaluation at render-time requires following back a number of times along the velocity field. Fine details that don’t diffuse out are one feature of the Gridless Advection approach that is useful in certain circumstances and isn’t too hard to capture approximately by making use of a construction called the Characteristic Map.

The mapping function phrases the movement of density in a special way. Instead of moving density around on a grid, we have a “mapping function” which changes at each time step that maps where material “came from” after a period of flow. The paper explains the details behind this map, but implementation follows directly from the fact that a map at time $n$ is based on the previous map:


$\B X_{n+1}(\B x) = \B X_n(\B x – \B u \Delta t)$


with the initial condition $\B X_0(\B x) = \B x$. With this map, the density at time $t$ is given by a combination of the map at time $t$ and the original density only:


$\rho_n(\B x) = \rho_0(\B X_n(\B x))$


In the Construct,  this looks like

ScalarField initial_density = ...;
VectorField X = identity();
for(...) {
  // Update the mapping function by advecting coordinates.
  // Please see the paper for details
  X = advect(X, velocity, dt);
  X = writeToGrid(X, identity(), dt);


  // Now the density is given by mapping the density through
  // the Characteristic Map
  auto render_density = warp(initial_density, X);
  render(render_density, ...);

A comparison of these techniques is shown below:


On the left is the vanilla traditional rendering of density, in the middle is the CM-based rendering of the same density, and on the right is the CM-based rendering of a different density (at no additional cost). As you can see, the CM-based approaches have much less diffuse appearance to the density.

 More to Come

More examples are coming! If you have questions or comments, please leave them here.

5 thoughts on “Extending the Simple Fluid Simulation”

  1. Hello Brandon. Thanx for next article and greate example.
    If you currently can not include code in Construct for volume rendering, then send me on email code for study, of course, if you can.
    Ok. How about a more complete example, where you can calculate not only the density, but also other information (temperature, energy and so on) and get fire with smoke and rendering data, using blackbody spectrum?
    And second. How about adding detail to the data obtained with low resolution and as a result to get – high resolution, without intensive computations (as an example of – wavelet turbulence for fluid, Theodore Kim). Maybe you have other studies in this area?

  2. Hi

    I was wondering if you have any tutorials using SPH 2D water simulation?
    Would be wonderful if there would be some clear explanations of how SPH works from a programmers point of view. Allt the variables being used, the calculations and so on.
    I’ve tried to find some good information or examples but all seem to be quite complicated with practically no explanations of why they do it the way the do.

    Any pointers to give out?

  3. Hi Brandon

    My comment seems to have been lost.
    I was wondering if you have any tutorial with source code how to simulate water in 2D?
    I’m trying to see how that would look in a platform game.

    Having a look at some solutions I seem to find that they calculate something like this.

    u = u / mass does that mean that they are trying to find the Unit Vector of U ?

Leave a Reply

Your email address will not be published. Required fields are marked *