Jay Taylor's notes

back to listing index

Fluid Simulation (with WebGL demo) - Zero Wind :: Jamie Wong

[web search]

Original source (jamie-wong.com)

Tags: math three.js webgl fluid-simulation computer-graphics navier-stokes fluid-dynamics jamie-wong.com

Clipped on: 2016-08-27

Fluid Simulation (with WebGL demo)

Posted on August 5, 2016

Click and drag to change the fluid flow. Double click to reset.

Note: The demos in this post rely on WebGL features that might not be implemented in mobile browsers.

About a year and a half ago, I had a passing interest in trying to figure out how to make a fluid simulation. At the time, it felt just a bit out of my reach, requiring knowledge of shaders, vector calculus, and numerical computation that were all just a little bit past my grasp. At the time, I was working through the Fluid Simulation Course Notes from SIGGRAPH 2007, and was struggling with the math. Now armed with a bit more knowledge and a lot more time, and with the help of other less dense resources like GPU Gems Chapter 38. Fast Fluid Dynamics Simulation on the GPU, I was finally able to figure out enough to get something working. I am still a beginner at simulations like this, and I’m going to brazenly ignore things like numerical stability, but hopefully I can help leapfrog you past a few places I got stuck.

We’re going to work with the simplest 2D fluid simulation, where the entire area is full of fluid, and we’re going to ignore viscosity.

The Velocity Field

As compared to a rigid, unrotating solid, where every bit of the thing has to be moving in the same direction at the same speed, each bit of a fluid might be moving differently. One way to model this is to use a vector field representing velocity. For any given (x,y) (x, y) (x,y) coordinate, this field will tell you the velocity of the fluid at that point.

u⃗(x,y)=(ux,uy) \vec u(x, y) = (u_x, u_y) u⃗(x,y)=(ux,uy)

A nice way to get an intuition about what a given field looks like is to sample the function in a grid of points, then draw arrows starting at each grid point whose size and orientation are dictated by the value of the function at that point. For the purposes of this post, we’re always going to be working over the domain x∈[−1,1] x \in [-1, 1] x∈[−1,1], y∈[−1,1] y \in [-1, 1] y∈[−1,1].

For instance, here’s a very simple field u⃗(x,y)=(1,0) \vec u(x, y) = (1, 0) u⃗(x,y)=(1,0) representing everything moving at a constant speed to the right.

And here’s a more interesting one u⃗(x,y)=(x,y) \vec u(x, y) = (x, y) u⃗(x,y)=(x,y) where things move away from the origin, increasing in speed the farther away from the origin they are.

We’re going to play with this one, u⃗(x,y)=(sin(2πy),sin(2πx)) \vec u(x, y) = \left( \sin (2 \pi y), \sin (2 \pi x) \right) u⃗(x,y)=(sin(2πy),sin(2πx)), since it creates some interesting visual results once we start making the fluid move accordingly.

For a more thorough introduction to vector fields, check out the Introduction to Vector Fields video on Khan Academy. The rest of the videos on multivariate calculus might prove useful for understanding concepts in fluid flow too.

Now then, let’s get things moving.

Advection

Advection is the transfer of a property from one place to another due to the motion of the fluid. If you’ve got some black dye in some water, and the water is moving to the right, then surprise surprise, the black dye moves right.

If the fluid is moving in a more complex manner, that black dye will get pulled through the liquid in a more complex manner.

Before we dive into how advection works, we need to talk a bit about the format of the data underlying these simulations.

The simulation consist of two fields: color and velocity. Each field is represented by a two dimensional grid. For simplicity, we use the same dimensions as the output pixel grid.

Previously, I described the velocity field as an analytical function u⃗(x,y) \vec u(x, y) u⃗(x,y). In practice, that analytical function is only used to initialize the grid values.

The simulation runs by stepping forward bit-by-bit in time, with the state of the color and velocity grids depending only on the state of the color and velocity grids from the previous time step. We’ll use u⃗(p⃗,t) \vec u(\vec p, t) u⃗(p⃗,t) to represent the velocity grid at 2d position p⃗ \vec p p⃗ and time t t t, and c⃗(p⃗,t) \vec c(\vec p, t) c⃗(p⃗,t) to represent the color in the same manner.

So how do we move forward in time? Let’s just talk about how the color field changes for now. If we consider each grid point as a little particle in the fluid, then one approach is to update the color of the fluid where that particle will be, one time step in the future.

c⃗(p⃗+u⃗(p⃗,t)Δt,t+Δt):=c⃗(p⃗,t) \vec c(\vec p + \vec u(\vec p, t) \Delta t, t + \Delta t) := \vec c(\vec p, t) c⃗(p⃗+u⃗(p⃗,t)Δt,t+Δt):=c⃗(p⃗,t)

In order to run these simulations in real-time at high resolution, we want to implement them on the GPU. It turns out that this method of updating the value at the new location of the particle is difficult to implement on the GPU.

First, the position we want to write, p⃗+u⃗(p⃗,t)Δt \vec p + \vec u(\vec p, t) \Delta t p⃗+u⃗(p⃗,t)Δt might not lie on a grid point, so we’d have to distribute the impact of the write across the surrounding grid points. Second, many of our imaginary particles might end up in the same place, meaning we need to analyze the entire grid before we decide what the new values of each grid point might be.

So, instead of figuring out where our imaginary particles at the grid points go to, we’ll figure out where they came from in order to calculate the next time step.

c⃗(p⃗,t+Δt):=c⃗(p⃗−u⃗(p⃗)Δt,t) \vec c(\vec p, t + \Delta t) := \vec c(\vec p - \vec u(\vec p) \Delta t, t) c⃗(p⃗,t+Δt):=c⃗(p⃗−u⃗(p⃗)Δt,t)

With this scheme, we only need to write to a single grid point, and we don’t need to consider the contributions of imaginary particles coming from multiple different places.

The last teensy hurdle is figuring out the value of c⃗(p⃗−u⃗(p⃗,t)Δt,t) \vec c(\vec p - \vec u(\vec p, t) \Delta t, t ) c⃗(p⃗−u⃗(p⃗,t)Δt,t), since p⃗−u⃗(p⃗,t) \vec p - \vec u(\vec p, t) p⃗−u⃗(p⃗,t) might not be at a grid point. We can hop this hurdle using bilinear interpolation on the surrounding 4 grid points (the ones linked by the dashed grey rectangle above).

Advecting the Velocity Field

Barring a bizarre sequence of perfectly aligned fans underneath the liquid, there’s no reason why the velocity field wouldn’t change over time. Just as black ink would move through the fluid, so too will the velocity field itself! Just as we can advect c⃗ \vec c c⃗ through u⃗ \vec u u⃗, we can also advect u⃗ \vec u u⃗ through itself!

Intuitively you can think of it this way: a particle moving in a certain direction will continue moving in that direction, even after it’s moved.

Since we’re storing velocity in a grid just like we did with color, we can use the exact same routine to advect velocity through itself. Below, watch the velocity change over time, with an initial velocity field of u⃗=(1,sin(2πy)) \vec u = (1, \sin(2 \pi y)) u⃗=(1,sin(2πy)).

See how the changes you make by dragging propagate through space via advection.

If you tried playing around with this, and saw a bunch of weird hard edges and might’ve thought to yourself “I don’t think fluids work like that…”, you’d be right. We’re missing an important ingredient, but before we look at the solution, let’s take a closer look at the problem.

Divergent Fields

Something about the velocity field below makes this intuitively not feel like a fluid. Fluids just don’t behave like this.

Same problem with this one…

If you look at where the arrows are pointing in each of the above 2 simulations, you’ll see that there are spots where the all the arrows point away from that spot, and others where all the arrows point toward that spot. Assuming the volume of the liquid is staying constant, the density of the fluid has to be changing for such a velocity field to be possible.

Water is roughly incompressible. That means that at every spot, you have to have the same amount of fluid entering that spot as leaving it.

Mathematically, we can represent this fact by saying a field is divergence-free. The divergence of a velocity field u⃗ \vec u u⃗, indicated with div(u⃗) div(\vec u) div(u⃗) or ∇⋅u⃗ \nabla \cdot \vec u ∇⋅u⃗, is a measure of how much net stuff is entering or leaving a given spot in the field. For our 2D velocity field, it’s defined like this:

∇⋅u⃗=(∂∂x,∂∂y)⋅(ux,uy)=∂ux∂x+∂uy∂y\begin{aligned} \nabla \cdot \vec u &= \left( \frac{\partial}{\partial x}, \frac{\partial}{\partial y} \right) \cdot \left( u_x, u_y \right) \\ &= \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \end{aligned}∇⋅u⃗=(∂x∂,∂y∂)⋅(ux,uy)=∂x∂ux+∂y∂uy

The first of the two not-very-fluidy fields above has an equation u⃗(x,y)=(x,y) \vec u(x, y) = (x, y) u⃗(x,y)=(x,y). Taking the divergence, we find:

∇⋅u⃗=∂∂x(x)+∂∂y(y)=1+1=2\begin{aligned} \nabla \cdot \vec u &= \frac{\partial}{\partial x}(x) + \frac{\partial}{\partial y}(y) \\ &= 1 + 1 \\ &= 2 \end{aligned}∇⋅u⃗=∂x∂(x)+∂y∂(y)=1+1=2

This positive value tells us that, in all places, more stuff is leaving that point than entering it. In physical terms, this means that the density is decreasing uniformly everywhere.

The other not-very-fluidy field has an equation u⃗(x,y)=sin(2πx),0) \vec u(x, y) = \sin(2 \pi x), 0) u⃗(x,y)=sin(2πx),0). If we look at its divergence, we see:

∇⋅u⃗=∂∂x(sin(2πx))+∂∂y(0)=2πcos(2πx)\begin{aligned} \nabla \cdot \vec u &= \frac{\partial}{\partial x}(\sin(2 \pi x)) + \frac{\partial}{\partial y}(0) \\ &= 2 \pi \cos (2 \pi x) \end{aligned}∇⋅u⃗=∂x∂(sin(2πx))+∂y∂(0)=2πcos(2πx)

Which tells us that in some places, density is increasing (where ∇⋅u⃗<0 \nabla \cdot \vec u < 0 ∇⋅u⃗<0), and in others, density is decreasing (where ∇⋅u⃗>0 \nabla \cdot \vec u > 0 ∇⋅u⃗>0).

Doing the same operation on the more fluidy looking swirly velocity field u⃗=(sin(2πy),sin(2πx) \vec u = (\sin ( 2 \pi y), \sin ( 2 \pi x ) u⃗=(sin(2πy),sin(2πx) that you saw in the section about advection, we discover ∇⋅u⃗=0 \nabla \cdot \vec u = 0 ∇⋅u⃗=0.

An incompressible fluid will have a divergence of zero everywhere. So, if we want our simulated fluid to look kind of like a real fluid, we better make sure it’s divergence-free.

Since our velocity field undergoes advection and can be influenced by clicking and dragging around the fluid, having an initially divergence-free velocity field isn’t enough to guarantee that the field will continue to be divergence-free. For example, if we take our swirly simulation and start advecting the velocity field through itself, we end up with something divergent:

So we need a way of taking a divergent field and making it divergence-free. To understand what force makes that happen in the real world, we need to talk about some honest-to-goodness physics.

Navier-Stokes

The Navier-Stokes equations describe the motion of fluids. Here are the Navier-Stokes equations for incompressible fluid flow:

∂u⃗∂t=−u⃗⋅∇u⃗−1ρ∇p+ν∇2u⃗+F⃗∇⋅u⃗=0\begin{aligned} & \frac{\partial \vec u}{\partial t} = -\vec u \cdot \nabla \vec u -\frac{1}{\rho} \nabla p + \nu \nabla^2 \vec u + \vec F \\ \\ & \nabla \cdot \vec u = 0 \end{aligned}∂t∂u⃗=−u⃗⋅∇u⃗−ρ1∇p+ν∇2u⃗+F⃗∇⋅u⃗=0

Where u⃗ \vec u u⃗ is the velocity field, ρ \rho ρ is density, p p p is pressure, ν \nu ν is the kinematic viscosity, and F⃗ \vec F F⃗ is external forces acting upon the fluid.

Since we’re pretending the viscosity of our fluid is zero, we can drop the ν \nu ν term in the first equation. In our simple simulation, external forces are only applied by dragging the mouse, so we’ll ignore that term for now, opting to allow it to influence the velocity field directly.

Dropping those terms, we’re left with the following:

∂u⃗∂t=−u⃗⋅∇u⃗−1ρ∇p \frac{\partial \vec u}{\partial t} = -\vec u \cdot \nabla \vec u -\frac{1}{\rho} \nabla p ∂t∂u⃗=−u⃗⋅∇u⃗−ρ1∇p

We can expand this to its partial derivative form, expanding vector components to leave us with only scalar variables.

[∂ux∂t∂uy∂t]=−[∂ux∂x∂ux∂y∂uy∂x∂uy∂y][uxuy]−1ρ[∂p∂x∂p∂y][∂ux∂t∂uy∂t]=[−ux∂ux∂x−uy∂ux∂y−1ρ∂p∂x−ux∂uy∂x−uy∂uy∂y−1ρ∂p∂y]\begin{aligned} \begin{bmatrix} \frac{\partial u_x}{\partial t} \\ \\ \frac{\partial u_y}{\partial t} \end{bmatrix} &= - \begin{bmatrix} \frac{\partial u_x}{\partial x} & \frac{\partial u_x}{\partial y} \\ \\ \frac{\partial u_y}{\partial x} & \frac{\partial u_y}{\partial y} \end{bmatrix} \begin{bmatrix} u_x \\ \\ u_y \end{bmatrix} - \frac{1}{\rho} \begin{bmatrix} \frac{\partial p}{\partial x} \\ \\ \frac{\partial p}{\partial y} \end{bmatrix} \\ \begin{bmatrix} \frac{\partial u_x}{\partial t} \\ \\ \frac{\partial u_y}{\partial t} \end{bmatrix} &= \begin{bmatrix} - u_x \frac{\partial u_x}{\partial x} - u_y \frac{\partial u_x}{\partial y} - \frac{1}{\rho} \frac{\partial p}{\partial x} \\ \\ - u_x \frac{\partial u_y}{\partial x} - u_y \frac{\partial u_y}{\partial y} - \frac{1}{\rho} \frac{\partial p}{\partial y} \end{bmatrix} \end{aligned}⎣⎢⎢⎢⎡∂t∂ux∂t∂uy⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡∂t∂ux∂t∂uy⎦⎥⎥⎥⎤=−⎣⎢⎢⎢⎡∂x∂ux∂x∂uy∂y∂ux∂y∂uy⎦⎥⎥⎥⎤⎣⎡uxuy⎦⎤−ρ1⎣⎢⎢⎢⎡∂x∂p∂y∂p⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡−ux∂x∂ux−uy∂y∂ux−ρ1∂x∂p−ux∂x∂uy−uy∂y∂uy−ρ1∂y∂p⎦⎥⎥⎥⎤

Remembering that these fields are all functions on (x,y,t) (x, y, t) (x,y,t), we can approximate the partial derivatives with finite differences. For instance, we can approximate the partial derivative of ux u_x ux with respect to t t t like so:

∂ux∂t≈ux(x,y,t+Δt)−ux(x,y,t)Δt \frac{\partial u_x}{\partial t} \approx \frac{u_x(x, y, t + \Delta t) - u_x(x, y, t)}{\Delta t} ∂t∂ux≈Δtux(x,y,t+Δt)−ux(x,y,t)

Because the procedure ends up being the same for both components, we’ll focus on only the x x x component here. Applying finite differences to all of the partial derivatives, we have this:

ux(x,y,t+Δt)−ux(x,y,t)Δt=−ux(x,y,t)ux(x+ϵ,y,t)−ux(x−ϵ,y,t)2ϵ−uy(x,y,t)ux(x,y+ϵ,t)−ux(x,y−ϵ,t)2ϵ−1ρp(x+ϵ,y,t)−p(x−ϵ,y,t)2ϵ \begin{aligned} \frac{u_x(x, y, t + \Delta t) - u_x(x, y, t)}{\Delta t} =& -u_x(x, y, t) \frac{u_x(x + \epsilon, y, t) - u_x(x - \epsilon, y, t)}{2 \epsilon} \\ & -u_y(x, y, t) \frac{u_x(x, y + \epsilon, t) - u_x(x, y - \epsilon, t)}{2 \epsilon} \\ & -\frac{1}{\rho} \frac{p(x + \epsilon, y, t) - p(x - \epsilon, y, t)}{2 \epsilon} \end{aligned} Δtux(x,y,t+Δt)−ux(x,y,t)=−ux(x,y,t)2ϵux(x+ϵ,y,t)−ux(x−ϵ,y,t)−uy(x,y,t)2ϵux(x,y+ϵ,t)−ux(x,y−ϵ,t)−ρ12ϵp(x+ϵ,y,t)−p(x−ϵ,y,t)

Ultimately what we want is u⃗(x,y,t+Δt) \vec u(x, y, t + \Delta t) u⃗(x,y,t+Δt), which will tell us, for a given point, what the velocity will be at the next time step. So let’s solve for that by rearranging the big long formula above:

ux(x,y,t+Δt)=ux(x,y,t)−ux(x,y,t)Δtux(x+ϵ,y,t)−ux(x−ϵ,y,t)2ϵ−uy(x,y,t)Δtux(x,y+ϵ,t)−ux(x,y−ϵ,t)2ϵ−1ρΔtp(x+ϵ,y,t)−p(x−ϵ,y,t)2ϵ \begin{aligned} u_x(x, y, t + \Delta t) = & u_x(x, y, t) \\ & - u_x(x, y, t) \Delta t \frac{u_x(x + \epsilon, y, t) - u_x(x - \epsilon, y, t)}{2 \epsilon} \\ & -u_y(x, y, t) \Delta t \frac{u_x(x, y + \epsilon, t) - u_x(x, y - \epsilon, t)}{2 \epsilon} \\ & -\frac{1}{\rho} \Delta t \frac{p(x + \epsilon, y, t) - p(x - \epsilon, y, t)}{2 \epsilon} \end{aligned} ux(x,y,t+Δt)=ux(x,y,t)−ux(x,y,t)Δt2ϵux(x+ϵ,y,t)−ux(x−ϵ,y,t)−uy(x,y,t)Δt2ϵux(x,y+ϵ,t)−ux(x,y−ϵ,t)−ρ1Δt2ϵp(x+ϵ,y,t)−p(x−ϵ,y,t)

If you look at the first three terms in this expression, what does it look like they conceptually represent? It looks like they represent the next velocity after we’ve taken into account changes due to the motion of the fluid itself. That sounds an awful like advection as discussed earlier. In fact, it will work quite well if we substitute the velocity field after it’s undergone advection. We’ll call the advected velocity field u⃗a \vec u ^ a u⃗a. So now we have:

ux(x,y,t+Δt)=uxa(x,y,t)−1ρΔtp(x+ϵ,y,t)−p(x−ϵ,y,t)2ϵ \begin{aligned} u_x(x, y, t + \Delta t) = u^a_x(x, y, t) -\frac{1}{\rho} \Delta t \frac{p(x + \epsilon, y, t) - p(x - \epsilon, y, t)}{2 \epsilon} \end{aligned} ux(x,y,t+Δt)=uxa(x,y,t)−ρ1Δt2ϵp(x+ϵ,y,t)−p(x−ϵ,y,t)

So after all of that, we have an equation that relates the velocity field at the next time tick to the current velocity field after it’s undergone advection, followed by application of pressure.

We know that a divergence-free field that undergoes advection isn’t necessarily still divergence-free, and yet we know that the Navier-Stokes equations for impressible flow represent divergence-free velocity fields, so therefore we have our answer about what in nature prevents the velocity field from becoming divergent: pressure!

Solving for Pressure

Now that we have an equation that relates u⃗ \vec u u⃗ to p p p. This is where the math gets messy. We start from the second Navier-Stokes equation for incompressible flow, applied at time t+Δt t + \Delta t t+Δt, and apply finite differences again:

∇⋅u⃗=0∂ux∂x+∂uy∂y=0ux(x+ϵ,y,t+Δt)−ux(x−ϵ,y,t+Δt)2ϵ+uy(x,y+ϵ,t+Δt)−uy(x,y−ϵ,t+Δt)2ϵ=0 \begin{aligned} \nabla \cdot \vec u & = 0 \\ \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} & = 0 \\ \frac{ u_x(x + \epsilon, y, t + \Delta t) - u_x(x - \epsilon, y, t + \Delta t) }{ 2 \epsilon } + \\ \\ \frac{ u_y(x, y + \epsilon, t + \Delta t) - u_y(x, y - \epsilon, t + \Delta t) }{ 2 \epsilon } & = 0 \end{aligned} ∇⋅u⃗∂x∂ux+∂y∂uy2ϵux(x+ϵ,y,t+Δt)−ux(x−ϵ,y,t+Δt)+2ϵuy(x,y+ϵ,t+Δt)−uy(x,y−ϵ,t+Δt)=0=0=0

Here, we can substitute our equations for u⃗ \vec u u⃗ expressed in terms of u⃗a \vec u ^ a u⃗a and p p p to get this monster:

0=12ϵ((uxa(x+ϵ,y,t)−1ρΔtp(x+2ϵ,y,t)−p(x,y,t)2ϵ)−(uxa(x−ϵ,y,t)−1ρΔtp(x,y,t)−p(x−2ϵ,y,t)2ϵ)+(uya(x,y+ϵ,t)−1ρΔtp(x,y+2ϵ,t)−p(x,y,t)2ϵ)−(uya(x,y−ϵ,t)−1ρΔtp(x,y,t)−p(x,y−2ϵ,t)2ϵ))\begin{aligned} 0 = \frac{1}{2 \epsilon} \left( \left( u^a_x(x + \epsilon, y, t) -\frac{1}{\rho} \Delta t \frac{p(x + 2\epsilon, y, t) - p(x, y, t)}{2 \epsilon} \right) \right. \\ \left. - \left( u^a_x(x - \epsilon, y, t) -\frac{1}{\rho} \Delta t \frac{p(x, y, t) - p(x - 2\epsilon, y, t)}{2 \epsilon} \right) \right. \\ \left. + \left( u^a_y(x, y + \epsilon, t) -\frac{1}{\rho} \Delta t \frac{p(x, y + 2\epsilon, t) - p(x, y, t)}{2 \epsilon} \right) \right. \\ \left. - \left( u^a_y(x, y - \epsilon, t) -\frac{1}{\rho} \Delta t \frac{p(x, y, t) - p(x, y - 2\epsilon, t)}{2 \epsilon} \right) \right) \end{aligned}0=2ϵ1((uxa(x+ϵ,y,t)−ρ1Δt2ϵp(x+2ϵ,y,t)−p(x,y,t))−(uxa(x−ϵ,y,t)−ρ1Δt2ϵp(x,y,t)−p(x−2ϵ,y,t))+(uya(x,y+ϵ,t)−ρ1Δt2ϵp(x,y+2ϵ,t)−p(x,y,t))−(uya(x,y−ϵ,t)−ρ1Δt2ϵp(x,y,t)−p(x,y−2ϵ,t)))

Rearranging to have all of the p p p terms on the left and all the u⃗a \vec u ^ a u⃗a terms on the right, and multiplying both sides by 2ϵ 2 \epsilon 2ϵ, we have:

−Δt2ϵρ(4p(x,y,t)−p(x+2ϵ,y,t)−p(x−2ϵ,y,t)−p(x,y+2ϵ,t)−p(x,y−2ϵ,t))=uxa(x+ϵ,y,t)−uxa(x−ϵ,y,t)+uya(x,y+ϵ,t)−uya(x,y−ϵ,t) -\frac{\Delta t}{2 \epsilon \rho} \left( \begin{matrix} 4 p(x, y, t) \\ -p(x + 2 \epsilon, y, t) \\ -p(x - 2 \epsilon, y, t) \\ -p(x, y + 2 \epsilon, t) \\ -p(x, y - 2 \epsilon, t) \end{matrix} \right) = \begin{matrix} u^a_x(x + \epsilon, y, t) \\ -u^a_x(x - \epsilon, y, t) \\ +u^a_y(x, y + \epsilon, t) \\ -u^a_y(x, y - \epsilon, t) \end{matrix} −2ϵρΔt⎝⎜⎜⎜⎜⎛4p(x,y,t)−p(x+2ϵ,y,t)−p(x−2ϵ,y,t)−p(x,y+2ϵ,t)−p(x,y−2ϵ,t)⎠⎟⎟⎟⎟⎞=uxa(x+ϵ,y,t)−uxa(x−ϵ,y,t)+uya(x,y+ϵ,t)−uya(x,y−ϵ,t)

Note: the above expression is a scalar expression, despite being laid out in a somewhat vector-y form.

At this point, it’s helpful to remember that, for the purposes of the simulation, we’re not interested in knowing the value of p p p everywhere: we only care about knowing its value at enough places to calculate the value of the velocities at the grid points. To meet that end, we can similarly calculate p p p on the grid. To accomplish this, we can make ϵ \epsilon ϵ the distance between adjacent grid cells.

The above equation yields a new equation for every (x,y) (x, y) (x,y) of a grid point we substitute. For the purposes of discussion, let’s assume that the gap between adjacent cells is 0.1 units, so ϵ=0.1 \epsilon = 0.1 ϵ=0.1. Let’s examine what the equation yields for (x,y)=(0.3,0.7) (x, y) = (0.3, 0.7) (x,y)=(0.3,0.7).

−Δt2ϵρ(4p(0.3,0.7,t)−p(0.3+2(0.1),0.7,t)−p(0.3−2(0.1),0.7,t)−p(0.3,0.7+2(0.1),t)−p(0.3,0.7−2(0.1),t))=uxa(0.3+0.1,0.7,t)−uxa(0.3−0.1,0.7,t)+uya(0.3,0.7+0.1,t)−uya(0.3,0.7−0.1,t)Δt2ϵρ(4p(0.3,0.7,t)−p(0.5,0.7,t)−p(0.1,0.7,t)−p(0.3,0.9,t)−p(0.3,0.5,t))=uxa(0.4,0.7,t)−uxa(0.2,0.7,t)+uya(0.3,0.8,t)−uya(0.3,0.6,t)\begin{aligned} -\frac{\Delta t}{2 \epsilon \rho} \left( \begin{matrix} 4 p(0.3, 0.7, t) \\ -p(0.3 + 2(0.1), 0.7, t) \\ -p(0.3 - 2(0.1), 0.7, t) \\ -p(0.3, 0.7 + 2(0.1), t) \\ -p(0.3, 0.7 - 2(0.1), t) \end{matrix} \right) & = \begin{matrix} u^a_x(0.3 + 0.1, 0.7, t) \\ -u^a_x(0.3 - 0.1, 0.7, t) \\ +u^a_y(0.3, 0.7 + 0.1, t) \\ -u^a_y(0.3, 0.7 - 0.1, t) \end{matrix} \\ \frac{\Delta t}{2 \epsilon \rho} \left( \begin{matrix} 4 p(0.3, 0.7, t) \\ -p(0.5, 0.7, t) \\ -p(0.1, 0.7, t) \\ -p(0.3, 0.9, t) \\ -p(0.3, 0.5, t) \end{matrix} \right) & = \begin{matrix} u^a_x(0.4, 0.7, t) \\ -u^a_x(0.2, 0.7, t) \\ +u^a_y(0.3, 0.8, t) \\ -u^a_y(0.3, 0.6, t) \end{matrix} \end{aligned}−2ϵρΔt⎝⎜⎜⎜⎜⎛4p(0.3,0.7,t)−p(0.3+2(0.1),0.7,t)−p(0.3−2(0.1),0.7,t)−p(0.3,0.7+2(0.1),t)−p(0.3,0.7−2(0.1),t)⎠⎟⎟⎟⎟⎞2ϵρΔt⎝⎜⎜⎜⎜⎛4p(0.3,0.7,t)−p(0.5,0.7,t)−p(0.1,0.7,t)−p(0.3,0.9,t)−p(0.3,0.5,t)⎠⎟⎟⎟⎟⎞=uxa(0.3+0.1,0.7,t)−uxa(0.3−0.1,0.7,t)+uya(0.3,0.7+0.1,t)−uya(0.3,0.7−0.1,t)=uxa(0.4,0.7,t)−uxa(0.2,0.7,t)+uya(0.3,0.8,t)−uya(0.3,0.6,t)

All the values on the right hand side of this equation are known, and on the left we have 5 unknowns: the value of p p p at 5 different grid locations.

If we repeat this process and evaluate (x,y) (x, y) (x,y) at every grid point, we get one equation with 5 unknowns for each grid location. If our grid has n×m n \times m n×m grid locations in it, then we have n×m n \times m n×m equations, each with 5 unknowns.

If you’re wondering about what’s happening at the edges, we’re going to lazily side-step that question by making our grid wrap around: if you ask for the velocity past the bottom edge, you’ll get a value near the top edge.

Before we move on, our notation is getting a bit clunky, so let’s clean it up a tad since we know we’re working on a grid. For the next part, we’ll say pi,j=p(iϵ,jϵ,t) p_{i,j} = p(i \epsilon, j \epsilon, t) pi,j=p(iϵ,jϵ,t), and we’ll stick all the known values together into a value d d d (for divergence), like so:

di,j=−2ϵρΔt(uxa((i+1)ϵ,jϵ,t)−uxa((i−1)ϵ,jϵ,t)+uya(iϵ,(j+1)ϵ,t)−uya(iϵ,(j−1)ϵ,t)) d_{i,j} = -\frac{2 \epsilon \rho}{\Delta t} \begin{pmatrix} u^a_x((i + 1) \epsilon, j \epsilon, t) \\ -u^a_x((i - 1) \epsilon, j \epsilon, t) \\ +u^a_y(i \epsilon, (j + 1) \epsilon, t) \\ -u^a_y(i \epsilon, (j - 1) \epsilon, t) \end{pmatrix} di,j=−Δt2ϵρ⎝⎜⎜⎛uxa((i+1)ϵ,jϵ,t)−uxa((i−1)ϵ,jϵ,t)+uya(iϵ,(j+1)ϵ,t)−uya(iϵ,(j−1)ϵ,t)⎠⎟⎟⎞

With this nicer notation, we can express the system of equations on pressure that we’re trying to solve like so:

4pi,j−pi+2,j−pi−2,j−pi,j+2−pi,j−2=di,j 4p_{i, j} - p_{i+2,j} - p_{i-2,j} - p_{i,j+2} - p_{i,j-2} = d_{i, j} 4pi,j−pi+2,j−pi−2,j−pi,j+2−pi,j−2=di,j

Iteratively Solving the Pressure Equation

Solving for pi,j p_{i, j} pi,j for every grid point analytically would be an enormous mess. Instead, we’re going to use an iterative method of solving this system of equations, where each iteration provides values closer and closer to a real solution. We’re going to use the Jacobi Method.

In the Jacobi method, we first rearrange our equation to isolate one term, like so:

pi,j=di,j+pi+2,j+pi−2,j+pi,j+2+pi,j−24 p_{i, j} = \frac{ d_{i, j} + p_{i+2,j} + p_{i-2,j} + p_{i,j+2} + p_{i,j-2} }{4} pi,j=4di,j+pi+2,j+pi−2,j+pi,j+2+pi,j−2

Next, we make an initial guess for all of our unknowns. We’ll call this initial guess pi,j(0) p_{i,j}^{(0)} pi,j(0), and just set it to 0 everywhere.

Here’s where the iteration comes in: our next guess, pi,j(1) p_{i, j}^{(1)} pi,j(1) is obtained by plugging in our initial guess into the above formula:

pi,j(1)=di,j+pi+2,j(0)+pi−2,j(0)+pi,j+2(0)+pi,j−2(0)4 p_{i, j}^{(1)} = \frac{ d_{i, j} + p_{i+2,j}^{(0)} + p_{i-2,j}^{(0)} + p_{i,j+2}^{(0)} + p_{i,j-2}^{(0)} }{4} pi,j(1)=4di,j+pi+2,j(0)+pi−2,j(0)+pi,j+2(0)+pi,j−2(0)

And, more generally, each iteration relies upon the previous one:

pi,j(k)=di,j+pi+2,j(k−1)+pi−2,j(k−1)+pi,j+2(k−1)+pi,j−2(k−1)4 p_{i, j}^{(k)} = \frac{ d_{i, j} + p_{i+2,j}^{(k-1)} + p_{i-2,j}^{(k-1)} + p_{i,j+2}^{(k-1)} + p_{i,j-2}^{(k-1)} }{4} pi,j(k)=4di,j+pi+2,j(k−1)+pi−2,j(k−1)+pi,j+2(k−1)+pi,j−2(k−1)

You would usually run this until the values of one iteration are equal to the values from the previous iteration, rounded to a certain accuracy. For our purposes, we’re more interested in this running in a consistent period of time, so we’ll arbitrarily run this for 10 iterations, and hope the result is accurate enough to look realistic.

For a bit of intuition on why this converges to a solution, check out Algebraic Pavel’s answer on Math Exchange.

All Together Now!

Phew! That was a lot to get through. Now let’s put it all together. Roughly, as pseudo-code, here’s our whole simulation:

initialize color field, c
initialize velocity field, u

while(true):
    u_a := advect field u through itself
    d := calculate divergence of u_a
    p := calculate pressure based on d, using jacobi iteration
    u := u_a - gradient of p
    c := advect field c through velocity field u
    draw c
    wait a bit

Here are the key formulas for those steps on grid coordinates (i,j) (i, j) (i,j), uncluttered by derivations:

Advecting field vecu vec u vecu through itself:

u⃗i,ja=u⃗a(x:=iϵ,y:=jϵ,t+Δt):=u⃗(x−ux(x,y)Δt,y−uy(x,y)Δt,t) \vec u^a_{i,j} = \vec u^a(x:=i \epsilon, y:= j \epsilon, t + \Delta t) := \vec u(x - u_x(x, y) \Delta t, y - u_y(x, y) \Delta t, t) u⃗i,ja=u⃗a(x:=iϵ,y:=jϵ,t+Δt):=u⃗(x−ux(x,y)Δt,y−uy(x,y)Δt,t)

Divergence of u⃗a \vec u_a u⃗a (multiplied by constant terms):

di,j=−2ϵρΔt(uxi+1,ja−uxi−1,ja+uyi,j+1a−uyi,j−1a) d_{i,j} = -\frac{2 \epsilon \rho}{\Delta t} ( u^a_{x_{i+1, j}} - u^a_{x_{i-1, j}} + u^a_{y_{i, j+1}} - u^a_{y_{i,j-1}} ) di,j=−Δt2ϵρ(uxi+1,ja−uxi−1,ja+uyi,j+1a−uyi,j−1a)

Pressure calculation Jacobi iteration step, with pi,j(0)=0 p_{i, j}^{(0)} = 0 pi,j(0)=0:

Subtracting the pressure gradient from the advected velocity field:

uxi,j:=uxi,ja−Δt2ρϵ(pi+1,j−pi−1,j)uyi,j:=uyi,ja−Δt2ρϵ(pi,j+1−pi,j−1) \begin{aligned} u_{x_{i, j}} &:= u^a_{x_{i,j}} -\frac{\Delta t}{2 \rho \epsilon} (p_{i + 1, j} - p_{i - 1, j}) \\ u_{y_{i, j}} &:= u^a_{y_{i, j}} -\frac{\Delta t}{2 \rho \epsilon} (p_{i, j + 1} - p_{i, j - 1}) \end{aligned} uxi,juyi,j:=uxi,ja−2ρϵΔt(pi+1,j−pi−1,j):=uyi,ja−2ρϵΔt(pi,j+1−pi,j−1)

Advecting the color field through the final velocity field:

ca⃗=ca⃗(x:=iϵ,y:=jϵ,t+Δt):=c⃗(x−ux(x,y)Δt,y−uy(x,y)Δt) \vec {c^a} = \vec{c^a}(x:=i \epsilon, y:= j \epsilon, t + \Delta t) := \vec c(x - u_x(x, y) \Delta t, y - u_y(x, y) \Delta t) ca⃗=ca⃗(x:=iϵ,y:=jϵ,t+Δt):=c⃗(x−ux(x,y)Δt,y−uy(x,y)Δt)

Implementation

Pulling all those steps together, you can make something like this! Woohoo! When I got this working for the first time, I was pretty ecstatic.

I won’t delve too far into the implementation, but you can have a look at it yourself: fluid-sim.js. It relies upon the elegant lightgl.js, which is an abstraction layer on top of WebGL that makes it much nicer to work with. Unlike THREE.js, it doesn’t make any assumptions about you wanting any concept of a camera or lighting or that you’re working in 3D at all.

The key technique for running the simulation efficiently is doing all the hard work on the GPU. To meet this need, all of the computations are done via the render to texture technique, ping-ponging which texture is being rendered to facilitate reading and writing to the same conceptual texture (e.g. reading from the velocity field and writing to the velocity field representing the next time step).

Each one of the major components of the algorithm is implemented in a separate shader. There’s a shader for advection, a shader for calculating the divergence, one for a single iteration of the Jacobi method, and another for subtracting the pressure gradient from the advected velocity.

References

To make this, I had to draw from a lot of difference references, many of which are linked inline in the post.

Fluid Simulation Course Notes from SIGGRAPH 2007: Now a textbook, this is a pretty mathematically dense tutorial. It took me 4 or 5 times reading through most sections to make sense of it, and ultimately I only understood parts of it after I did the derivations myself. It uses the more complex conjugate gradient method instead of the Jacobi method to solve the system of pressure equations, which I got completely lost in, and abandoned. It delves into a lot of arguments about numerical accuracy and uses a more complex grid layout than I did, which I still don’t follow fully. It also has resources for other kinds of fluid simulations, like heightfield simulation, and smoothed particle hydrodynamics.
GPU Gems Chapter 38. Fast Fluid Dynamics Simulation on the GPU: This was the single most useful reference I found, and describes something very similar to this post. It walks through specific implementation ideas, and gave me a much better intuition for advection. Some of the math (or at least the notation) seems shaky here. I think the Gaussian “splat” formula is missing a negative sign inside of the exp() exp() exp(), and I’m not sure what the notation (u⃗⋅∇)ux (\vec u \cdot \nabla) u_x (u⃗⋅∇)ux means in the first Navier-Stokes equation, since ∇⋅u⃗=0 \nabla \cdot \vec u = 0 ∇⋅u⃗=0 in the second equation.
“Elementary Linear Algebra” by Ron Larson, Section 10.2: Iterative Methods for Solving Linear Systems. This had a much clearer explanation of the Jacobi method than the GPU Gems chapter that allowed me to derive the pressure solve iteration myself. The full textbook can be found here: “Elementary Linear Algebra” on amazon.com.
Jonas Wagner’s fluid simulation on canvas, and particularly the source for it (fluid.js) were helpful for understanding what a full solution actually looks like. It’s also how I found the GPU Gems article in the first place. Jonas went on later to reimplement his solution in WebGL: WebGL Fluid Simulation.