Simulating Fluids, Fire, and Smoke in Real-Time
Notes on the math, algorithms, and methods involved in simulating fluids like fire and smoke in real-time.
Source code for this article can be found on my GitHub.
Fire is an interesting graphics problem. Past approaches generally faked it. For example,
Lord of the Rings
But in the last 10 years GPUs have made fast fluid simulation easy.
Basic fluid dynamics algorithms are straightforward to implement on the GPU
This post takes notes on how fire can be simulated on the GPU. It walks through the math behind fluid dynamics, parallel algorithms for modeling fluids, and the extra combustion bits that make fire special. Readers should have a reasonable background in vector calculus and differential equations (know how to take the gradient of a vector). Demos are implemented with WebGL.
1. Fluid Simulation
Before we simulate fire, we need to simulate fluid. We assume our fluid is incompressible and inviscid, which will vastly simplify our problem.
1.1 Basic Fluid Dynamics
Suppose \(D\) is a region in space filled with a fluid. At any point \( \mathbf{x} \in D \) and time \(t\), the fluid has velocity \(\mathbf{u}(\mathbf{x}, t)\). Computationally, we can represent the 2D velocity field \( \mathbf{u} \) with an \( N \times N \) grid, where the equally spaced grid points give the value of the velocity field at that point in space.
Ex: A 16\(\times\)16 grid representing \( \mathbf{u} = (y, -x) \)
What will happen if we put a drop of dye in the fluid?
Let's define a scalar field \( \psi (\mathbf{x}, t) \) as the density of the dye at any point in space and time. The transport of quantities like \( \psi \) within a fluid by the fluid's velocity is called advection. Given some fluid's velocity field and an initial density field of our dye, we'd like to see how the dye's density everywhere evolves over time by simulating its advection through the fluid.
A Naive Method for Advection
One idea
But this is tricky to parallelize, since 2 grid points can end up in the same target point after forward evaluation. And in practice, the target point will fall between grid points, meaning it has to be interpolated into the surrounding grid points. Finally, this process is unstable for time steps above some number, causing \( \psi \) to blow up.
The Advection Partial Differential Equations
This whole time we've been solving a partial differential equation! If we're going to derive a stable method for advection, we'll need to first get an explicit expression for this PDE. Let's start from first principles.
Consider a fixed region of space \(W\) (that is, \(W\) does not vary with time). The total mass of dye within \(W\) is \( \int_{W} \psi dV\). Over time, the change in mass is: $$ \frac{d}{dt} \int_{W} \psi (\mathbf{x}, t) dV = \int_{W} \frac{\partial}{\partial t} \psi (\mathbf{x}, t) dV \\ $$ Now, letting \(S\) denote the surface of \(W\) and \( \mathbf{n} \) the outward normal vector defined along the surface, we can examine the mass flow rate through the surface of \(W\). In particular, observe that the volume flow rate - the volume of fluid that flows through per second - across \(S\) per unit area is \( \mathbf{u} \cdot \mathbf{n} \) and the mass flow rate per unit area is \( \psi \mathbf{u} \cdot \mathbf{n} \).
This gives us the law of conservation of mass in integral form: $$ \frac{d}{dt} \int_{W} \psi dV = - \int_{S} \psi \mathbf{u} \cdot \mathbf{n} dA $$ Can we get rid of the integrals and say something similar for points? By divergence theorem, the above is equivalent to $$ \int_{W} [\frac{\partial \psi}{\partial t} + \nabla \cdot (\psi \mathbf{u})] dV = 0 $$ Then for a unit subregion \( W = dV \), we can say that $$ \frac{\partial \psi}{\partial t} + \nabla \cdot (\psi \mathbf{u}) = 0 $$ This gives us an explicit PDE that we need to solve for \( \psi \)!
Hmm... we could stop here, but we might be able to simplify this more. Specifically, it looks like we could isolate out a term \( \nabla \cdot \mathbf{u} \) that goes to zero because of incompressibility. $$ \begin{aligned} &\frac{\partial \psi}{\partial t} = - \nabla \cdot (\psi \mathbf{u}) \\ &= - (\frac{\partial}{\partial x} \psi u + \frac{\partial}{\partial y} \psi v) \\ &= - (\frac{\partial \psi}{\partial x} u + \frac{\partial u}{\partial x} \psi + \frac{\partial \psi}{\partial y} v + \frac{\partial v}{\partial y} \psi) \\ &= - (\psi \nabla \cdot \mathbf{u} + \mathbf{u} \cdot \nabla \psi) \\ &= - \mathbf{u} \cdot \nabla \psi \\ \end{aligned} $$ Applying our incompressibility constraint \( \nabla \cdot \mathbf{u} = 0 \) at the end yields a scalar PDE, the first of our incompressible flow advection equations: $$ \frac{\partial \psi}{\partial t} = \text{advection} (\mathbf{u}, \psi) = - \mathbf{u} \cdot \nabla \psi \tag{1} $$ $$ \frac{\partial \mathbf{v}}{\partial t} = \text{advection} (\mathbf{u}, \mathbf{v}) = - \mathbf{u} \cdot \nabla \mathbf{v} \tag{2} $$ Eq. 2 for advecting a vector field \( \mathbf{v} \) through our velocity field can be derived similarly to the scalar advection equation.
A Stable Method for Advection
Let's look closely at eqn. (1): $$ \frac{\partial \psi}{\partial t} = - \mathbf{u} \cdot \nabla \psi $$
Notice that the right-hand term is a directional derivative in the \( -\mathbf{u} \) direction. This gives us a wonderful new method for advecting \( \psi \) by an incompressible fluid - starting at a grid point \( \mathbf{x} \), trace the fluid velocity backwards, replacing the value at our original point by the value that we land on (if we land between points, we can interpolate): $$ \psi (\mathbf{x}, t + \Delta t) = \psi (\mathbf{x} - \mathbf{u} (\mathbf{x}, t), t) $$ In GPU pseudocode:
global Vec2Field u;
global FloatField density;
global float dt;
// Run for each point in our scalar grid that we want to update
float advectPoint(vec2 x) {
vec2 coord = x - dt * getVec2At(u, x);
return getFloatAt(density, coord);
}
This method is called
Semi-Lagrangian advection and was invented in 1999 by Jos Stam
- It's extremely easy to parallelize because each grid point only gets updated once per iteration.
- It's unconditionally stable. Why? Observe that for any grid point, the maximum value it can get updated to is the maximum value of all the grid points.
For a fixed velocity field fulfilling the incompressibility constraint, it works great.
Click anywhere above to drop some dye in the flow
1.2 The Navier-Stokes Equations
So far we've found a model that describes how scalar properties of a fluid evolve over time, assuming the flow is fixed. What about the fluid flow itself - how does the velocity field \( \mathbf{u} \) affect itself over time?
The Navier-Stokes equations
let u = createVectorGrid();
let density = createScalarGrid();
let p = createScalarGrid();
while (true) {
// Solve for the next velocity field.
u = advect(u, u);
p = computePressure(...);
u = u - gradient(p);
// Advect dye through the new velocity field.
density = advect(u, density);
}
Let's take a closer look at each of these.
Self-Advection
From our incompressible advection equations, we can see that the self-advection term is the advection of the fluid's velocity field \( \mathbf{u} \) through itself: $$ \text{advection} (\mathbf{u}, \mathbf{u}) = - \mathbf{u} \cdot \nabla \mathbf{u} \tag{4} $$
Where do the other terms come from? Well, advecting \( \mathbf{u} \) through itself yields a new velocity field \( \mathbf{u}^\prime \) (computable via the Semi-Lagrangian backtracing algorithm from above): $$ \mathbf{u}^\prime = \mathbf{u} - \mathbf{u} \cdot \nabla \mathbf{u} $$
Pressure
We don't know if this new velocity field follows the incompressibility constraint (e.g. has zero divergence). So the pressure term \( p \) needs to correct this somehow: $$ \nabla \cdot (\mathbf{u}^\prime - \nabla p) = 0 $$ We rearrange this to get $$ \nabla^2 p = \nabla \cdot \mathbf{u}^\prime \tag{5} $$ which is a type of equation known as a Poisson equation, where the left-hand side is the Laplacian of an unknown scalar field and the right-hand side is a known scalar. Solving this Poisson equation is really the slowest computational step in fluid simulation, for reasons we will see shortly.
Solving for Pressure
So how do we solve this particular PDE for \( p \)? Well, we know the value of our candidate velocity field \( \mathbf{u}^\prime \)
at all of our grid points, so we can approximately compute the right-hand side of the Poisson equation by applying a discrete
version
Then we can use a discrete version of the Laplacian $$ \nabla^2 p \approx p_{i+1, j} + p_{i-1, j} + p_{i, j+1} + p_{i, j-1} - 4p_{i, j} $$ to transform the whole equation into a linear equation with five unknowns.
But really, we are solving the Poisson equation (5) over all of space at once, so for an \( N \times N \) grid, we end up with a system of \( N^2 \) linear equations with exactly \( N^2 \) unknowns! So we end up with the familiar old equation $$ \mathbf{Ax} = \mathbf{b} $$ where \( \mathbf{A} \) is a matrix applying the Laplacian operator to the whole grid and \( \mathbf{b} \) is a vector containing the velocity field's divergence at all grid points.
There are many off-the-shelf algorithms for solving linear systems exactly. Unfortunately for us, even the fastest algorithms scale superlinearly with our grid size.
Solving for Pressure... Efficiently
If we're going to make a real-time simulation, we need to go fast. Can we leverage the GPU?
Well, it's not really possible to achieve an exact solution to the linear system efficiently, but we should note that the linear system is already an approximation to the Poisson equation. And it is possible to achieve arbitrarily accurate approximations with iterative methods - which begin with an estimate and improve solution accuracy every iteration - so we can just pick an iterative algorithm and run it until we have something that's "good enough". One particularly simple and easy-to-implement iterative algorithm for solving linear equations is the Jacobi method.
We start with the very first equation in the system: $$ A_{11}x_1 + A_{12}x_2 + ... + A_{1n}x_n = b_1 $$ At the \(k\)th iteration, given some guess \( \mathbf{x}^k \) for the solution \( \mathbf{x} \), we have some error. We can use this error to update our guess for \( x_1 \) as follows: $$ x_1^{k+1} = \frac{ b_1 - A_{12}x_2^k - ... - A_{1n}x_n^k }{ A_{11} } $$ In Jacobi, our guesses for all elements of \( \mathbf{x} \) are executed in parallel, giving a perfect match for implementation on the GPU. In pseudocode:
global FloatField divergence;
global FloatField pressure;
global float texelSize;
// Run for each point in our pressure grid that we want to update
float iterateJacobi(vec2 x) {
float div = getFloatAt(divergence, x);
float L = getFloatAt(pressure, x + vec2(-texelSize, 0.));
float R = getFloatAt(pressure, x + vec2(texelSize, 0.));
float T = getFloatAt(pressure, x + vec2(0., texelSize));
float B = getFloatAt(pressure, x + vec2(0., -texelSize));
return (div - L - R - T - B) / -4.;
}
It's worth noting that other, faster-converging solvers can also be
implemented on the GPU, like the Conjugate Gradient method and the Multigrid method. But depending on
the fluid and application, pressure accuracy may not be as important as advection accuracy or ease of implementation.
For smoke and fire, changes in fluid volume aren't as apparent as they are for fluids like water
Summary: Simulating Navier Stokes
The math behind Navier-Stokes can be a little bit dense, but at a high-level, simulating a fluid by solving the equations comes down to a few key update procedures on a grid per timestep. For our dye problem, here's our simulation might look:
let u = createVectorGrid();
let density = createScalarGrid();
let div = createScalarGrid();
let p = createScalarGrid();
while (true) {
// Solve for the next velocity field.
u = advect(u, u);
// Enforce incompressibility with pressure projection.
div = divergence(u);
for (let i = 0; i < JACOBI_ITERATIONS; i++) {
p = updatePressure(p, div);
}
u = u - gradient(p);
// Advect dye through the new velocity field.
density = advect(u, density);
}
1.3 Vorticity Confinement
Using a grid to store our velocity field is extremely convenient, but it results in unwanted numerical smoothing whenever we have to interpolate values between grid points. This combined with the relatively coarse approximation of a first-order Semi-Lagrangian advection scheme has the effect of dissipating out turbulent vortices in our flow. Physically, the velocity field loses energy, and the end result is generally overly smooth, "boring" fluid flow.
One way to combat lost vorticity is to increase the resolution of our grid, but this isn't really feasible
for real-time simulations that have limited computational resources. What we would ideally
like to do is find all the small details that get smoothed over each step of the simulation, and amplify
them. This process is called vorticity confinement - admittedly, it's not totally realistic, but succeeds in
preserving small scale details in more or less physically correct locations
The smallest turbulent features we can find are the vortices centered at each grid point in our simulation. We can measure the intensity of these vortices (the vorticity of them) by taking the curl of \( \mathbf{u} \) at each point, and amplify them by essentially adding a circular flow scaled by vorticity about each point. Mathematically, the vorticity is defined by $$ \bm{\omega} = \nabla \times \mathbf{u} $$ For each grid point, we compute a normalized location vector that points to the highest nearby vorticity concentration: $$ \mathbf{N} = \frac{ \nabla | \bm{\omega} | }{ | \nabla | \bm{\omega} | | } $$ And finally, we compute the confined vorticity vector field and add it to our flow: $$ \mathbf{f_{conf}} = \epsilon (\mathbf{N} \times \bm{\omega}) $$ $$ \mathbf{u_{conf}} = \mathbf{u} + \mathbf{f_{conf}} $$ Here, the confinement constant \( \epsilon > 0 \) is a parameter controlling the amount of small scale detail added back to the flow. Even low confinement levels (around 0-15) can make a huge difference, especially for simulations using Semi-Lagrangian advection schemes, and higher confinement levels can create highly stylized, billowing flows.
Click and drag to drop some dye in the turbulent simulation above
On the GPU, we can compute curl and confinement like so:
global Vec2Field u;
global float texelSize;
// Run to get curl for each point in grid
float computeCurl(vec2 x) {
float L = getVec2At(u, x + vec2(-texelSize, 0.)).y;
float R = getVec2At(u, x + vec2(texelSize, 0.)).y;
float T = getVec2At(u, x + vec2(0., texelSize)).x;
float B = getVec2At(u, x + vec2(0., -texelSize)).x;
return (R - L) - (T - B);
}
global Vec2Field curl;
global float confinement;
// Run to get confinement force for each point in grid
vec2 confinementForce(vec2 x) {
float L = getFloatAt(curl, x + vec2(-texelSize, 0.));
float R = getFloatAt(curl, x + vec2(texelSize, 0.));
float T = getFloatAt(curl, x + vec2(0., texelSize));
float B = getFloatAt(curl, x + vec2(0., -texelSize));
float C = getFloatAt(curl, x);
vec2 N = vec2(abs(T) - abs(B), abs(R) - abs(L));
N = N / length(N);
return confinement * C;
}
The full simulation with turbulence:
let u = createVectorGrid();
let density = createScalarGrid();
let div = createScalarGrid();
let p = createScalarGrid();
let curl = createVectorGrid();
while (true) {
// Solve for the next velocity field.
u = advect(u, u);
// Use vorticity confinement to amplify turbulence of velocity field.
curl = computeCurl(u);
u = u + confinementForce(curl, CONFINEMENT);
// Enforce incompressibility with pressure projection.
div = divergence(u);
for (let i = 0; i < JACOBI_ITERATIONS; i++) {
p = updatePressure(p, div);
}
u = u - gradient(p);
// Advect dye through the new velocity field.
density = advect(u, density);
}
Curl-Noise Turbulence
Curl noise is a method that essentially does the same thing as vorticity confinement, but instead of measuring and amplifying the vorticity of the velocity field, a scalar vorticity field is made from scratch using noise functions. Mathematically, we can combine vorticity confinement and curl-noise turbulence by synthesizing a random vorticity field $$ \bm{\phi} = \text{rand} * \mathbf{z} $$ Then computing our final vorticity field by $$ \bm{\omega}^* = \bm{\omega} + \bm{\phi} $$ Fast-moving, highly-turbulent fluids like smoke and fire benefit the most from vorticity confinement and curl noise, and in practice the curl noise field \( \bm{\phi} \) both evolves with time and is also advected by the fluid flow.
2. Fire Simulation
If you've gotten this far, pat yourself on the back! The methods in the previous section let us efficiently
and accurately simulate fluids with varying physical parameters (oil, water, honey)
assuming the fluid domain is a fixed space. Those interested in handling varying domains (that is, fluids that occupy different regions
within the grid, like a half-full cup of water that sloshes around) will want to explore accounting for different boundary conditions
within the grid simulation
Simulating fire and smoke requires a couple additions. First, we'll need to add channels representing fuel and temperature to our simulation, and model the combustion of fuel to create heat. Next we'll address how hot pockets of our fluid rise with a thermal buoyancy model, and finally, we'll need to render our flames correctly, taking into account blackbody radiation of the flames, human perception of light, and fire movement.
2.1 A Basic Combustion Model
Chemically, fire is caused by the oxidation of a fuel material in a reaction that releases both heat and light. In our case, we can assume that any fuel in our system has already ignited and is actively adding heat; we won't worry about the problem of unignited fuel.
To be more specific, let's define a scalar field \( \rho \) where \( 0 \leq \rho \leq 1 \) represents the density of fuel and another scalar field \( T > 0 \) representing the temperature of the fluid everywhere. At every timestep, temperature is added to the system by the fuel, which burns at a given burn temperature: $$ T^\prime = \text{max} ( T, \rho * T_{\text{burn}} ) $$ Of course, temperature isn't static - heat diffuses from hot to cold areas, and with fluids in particular, large-scale movements of molecules transport heat. The combination of these 2 processes defines heat convection, and conveniently, we already have a mathematical model for how it works - advection! Simulation-wise, we advect our temperature field along our velocity field. Since any reacting molecules are also moved by the fluid, we should advect fuel as well. The heat itself also affects the movement of the fluid - we'll see how to handle this shortly.
Furthermore, hot molecules radiate off temperature as light
To complete our combustion model, note that our fuel is always burning (we can imagine it as the density of ionized gas particles that give off thermal energy and return to a lower energy state), so every timestep we dissipate it by some given burn rate \( \gamma_{fuel} \): $$ \rho^\prime = \rho (1 - \gamma_{fuel})^{\Delta t} $$
2.2 Thermal Buoyancy
So far, our temperature field doesn't do anything to our fluid flow. But it should - hot pockets of air should expand and rise, and cooler pockets should fall. We can model this with a thermal buoyancy force. Since we're assuming incompressibility, we won't actually handle air expansion, but the fluid flow should experience an upward force depending on temperature: $$ \mathbf{u}^\prime = \mathbf{u} + (\beta T \Delta t) \mathbf{j} $$ Here, \( \beta \) is a given positive buoyancy constant, and \( \mathbf{j} \) is the upward unit vector.
Adding a combustion model and thermal buoyancy force gives us a fantastic simulator for a decidedly "fire-like" fluid - with the right values of buoyancy and cooling, we can get bulky, billowing plumes of material. Not exactly flames, but very similar to smoke.
Tap and drag in the simulation below to inject some combusting fuel. The displayed pixels represent density of smoke particles, which dissipate at a constant rate instead of being used up during combustion, but are still advected by the fluid simulation.
Click and drag to add smoke above
The simulation code builds off the basic fluid routines:
let u = createVectorGrid();
let density = createScalarGrid();
let div = createScalarGrid();
let p = createScalarGrid();
let curl = createVectorGrid();
let fuel = createScalarGrid();
let temp = createScalarGrid();
while (true) {
// Solve for the next velocity field.
u = advect(u, u);
// Combustion step.
temp = combust(temp, fuel);
// Use vorticity confinement to amplify turbulence of velocity field.
curl = computeCurl(u);
u = u + confineVorticity(curl, CONFINEMENT);
// Add thermal buoyancy.
u = u + buoyancy(temp);
// Enforce incompressibility with pressure projection.
div = divergence(u);
for (let i = 0; i < JACOBI_ITERATIONS; i++) {
p = updatePressure(p, div);
}
u = u - gradient(p);
// Advect density, fuel, and temperature through the new velocity field.
density = advect(u, density);
temp = advect(u, temp);
fuel = advect(u, fuel);
}
2.3 Fire Rendering
Fire is a participating medium,
meaning it emits light through blackbody radiation
The temperature-to-color spectrum as described by Planck's law
Planck's Law describes the spectral density of light radiated by a black body at a given temperature \( T \): $$ M(\lambda, T) = \frac{ c_1 }{ \lambda^5 } \frac{ 1 }{ \exp{ \frac{c_2}{\lambda T} } - 1 } $$ where $$ c_1 = 2 \pi h c^2 \\ c_2 = \frac{hc}{k} $$ and \(h\), \(c\), and \(k\) are Planck's constant, the speed of light, and Boltzmann's constant, respectively.
After implementing blackbody rendering using fragment shaders, we have a complete fire simulation!
Click and drag to add fire above
That's it for these notes! There is of course much more to fluid and fire simulation not covered here, like different (e.g. non-grid-based) techniques for solving the same problem of simulation within a fixed volume, different problems to solve involving varying domains or dynamic obstacles, enhancements to rendering like more accurate blackbody radiation, light scattering, or post-processing effects. Helpful introductions to these topics can be found in references below.