Simulation

How do we simulate Fire?

A common approach to the problem is to model the fire as a fluid with a buoyancy force greater than the gravitational force. While researching the approaches for fire simulation, I came across two approaches: particle based simulations and grid based simulations. I was initially drawn to a particle based simulation, specifically the position based fluids method developed by Macklin and Müller, but this approach would leave me to learn more advanced rendering techniques to transform particles to have fire characteristics, which wasn't a priority of mine. The alternative, grid based simulations, had more complexity for the simulation techniques, but the rendering techniques were more straightforward and more familiar to me. Additionally, grid based simulations have been researched and implemented more.

Grid based simulations model a velocity vector field, $\mathbf{u}(\mathbf{x}, t)$, at discrete grid locations. Each grid cell contains a vector representing the direction and speed of the fluid at that point in the simulation space.

Other scalar fields exist in our simulation like fuel, temperature, and smoke. These scalar fields are advected by the velocity field at each timestep of the simulation. Not only do these scalar quantities move through the vector field, but the vector field also advects itself. This gives the fluid simulation the realistic look we're looking for, but adds more complexity that we will tackle next.

The Navier Stokes Equations

The Navier Stokes Equations for incompressable flow describe how a fluid's vector field evolves over time at any given point.

Although very intimidating, this equation has been vastly simplified given our incompressibility constraint, $\nabla \cdot \mathbf{u} = 0$, and we can connect each term to a physical concept.

Advection

As explained earlier, advection is defined as the velocity of the fluid moving objects, densities, and other values along the field. This term represents the velocity field advecting itself.

Pressure

Forces don't instantly propagate through a fluid. Each molecule close to the source of the force pushes on other molecules adjacent to them, and those particles push on other molecules adjacent to them, and so on. This is pressure buildup and adds acceleration to the liquid.

Diffusion

There are different types of liquid, specifically those that vary in viscosity - honey moves much slower compared to water. This term represents how viscosity, or the ability for a fluid to resist flow, manipulates the velocity field.

External Forces

This term represents any forces outside the fluid that affect its velocity field.

Assumptions and Constraints

As mentioned earlier, incompressibility, $\nabla \cdot \mathbf{u} = 0$, is applied to the simulation. Additionally, this simulation will be inviscid, $\nu = 0$, so we don't have to consider the diffusion term. This is a beneficial trade off as fire doesn't visually behave like a viscous fluid.

Semi-Lagrangian Advection

Relative to other concepts in fluid simulation, advection seems relatively straightforward, but it has some hidden nuances. A naive approach would be to apply Euler's method - just move whatever quantity at position given by $\mathbf{r}$ with input of time $t$ by the magnitude of the velocity at time $t$ scaled by the timestep.

This approach is only stable in simulations with small timesteps, meaning the simulation won't "blow up". Having a small timestamp is computationally expensive, which limits doing the simulation in real time. This presents another issue because we are utilizing the GPU with parallelization. A cell in our simulation grid will inevitably have two or more quantities converge on it, creating data races.

The solution is to use a Semi-Lagrangian advection, an implicit method first introduced by Jos Stam in Stable Fluids, 1999 [1]. Instead, we backtrace where a quantity would be given the velocity, and move that quantity to the current position. To update a quantity $q$ for each timestep, we follow:

Where $q(\mathbf{x}, t)$ is the value of the quantity $q$ at position $\mathbf{x}$ at time $t$.

This solution is stable with any timestep and is easily parallelizable on the GPU. Although, this solution adds dissipation to the simulation, but this will be readjusted with vorticity confinement later.

Pressure Calculation

Now that we have a discrete and parallelized method of advection, we now need to apply the external forces and pressure terms to the vector field. Applying forces is straightforward, so I don't feel the need to talk about it, but pressure presents another story.

After applying advection and the external forces to our vector field for a discrete timestamp, we are left with a new vector field $u'$. This field isn't guaranteed to be incompressible (have zero divergence). Luckily, the Helmholtz-Hodge decomposition will be helpful in transforming a non-zero divergence vector field to a divergence-free vector field.

Let's define projection $\mathbb{P}$ that projects a vector field $u'$ onto it's divergence-free vector field $u$. If we apply that to the equation above, $\mathbb{P}\mathbf{u'} = \mathbb{P}\mathbf{u} + \mathbb{P}(\nabla p)$, and since by definition of $\mathbb{P}$, $\mathbb{P}\mathbf{u'} = \mathbb{P}\mathbf{u}$, we get that $\mathbb{P} (\nabla p) = 0$. Let's use these ideas with the Navier Stokes equation by applying the projection $\mathbb{P}$ to it.

Because $u$ is divergence-free, so is its derivative on the left hand side, so $\mathbb{P}\frac{\partial \mathbf{u}}{\partial t} = \frac{\partial \mathbf{u}}{\partial t}$. Additionally, $\mathbb{P}(\nabla p) = 0$, so we're left with:

This is a great insight because to solve the equation for each time stamp, we need to self-advect the velocity field, apply diffusion (in our case we are skipping this), apply external forces, and apply the projection $\mathbb{P}$ to ensure that the next iteration is divergence-free. To apply the projection in practice we can reuse the equation from Helmholtz-Hodge decomposition again by applying the divergence operator to both sides:

Since $\nabla \cdot u = 0$. We have a Poisson equation that we can use to solve for pressure. A key insight to note:

Applying the acceleration from the pressure term is exactly what is keeping the velocity field divergence free. Using the Helmholtz-Hodge decomposition theorem, we have discovered how to calculate pressure from the divergent velocity field and use it to create a divergence-free velocity field [4].

Solving Poisson Equations

All that is left now is to solve for pressure using the equation we defined above:

Since we are working with our discrete 3D grid, we need to use the finite difference form for the laplacian and divergence operator.

This assumes that the difference between each grid is unit length (h = 1). Now we plug them in:

We now have a linear equation with seven unknowns, but we want to solve for pressure at each of our $N^3$ grid locations. That leaves us with a system $N^3$ linear equations with $N^3$ unknowns because we want to solve for $i, j, k \in N^3$. It would take a lot of time to solve this system of linear equations exactly, so we can use the Jacobi method (which is well parallelized on the GPU) with $k=20$ iterations to find an approximate solution.

We now fully have a way to calculate the change in our velocity field for a discrete time step in our simulation!

Vorticity Confinement

Although we now have a stable simulation, we had to give up some things to achieve it. Smoke, air, or any gas you can think of has circular flows that spiral like pinwheels, which is referred to as vorticity. Fedkiw et al. explained how numerical dissipation caused by simulation on a coarse grid damps out these interesting features [2]. The process of adding these visualizations back is referred to as vorticity confinement. To do so, we first calculate the vorticity of our field, $\omega = \nabla \times u$, and use that to form a normalized vorticity vector field:

Where ${\eta} = \nabla | \omega |$. The resulting vector field, $\Psi$, consists of vectors that point towards the highest nearby vorticity concentration. We then use $\Psi$ with a confinement constant, $\epsilon$, to create an artificial force.

This force is then added to our velocity field to bring back the billowing effect of our fire and smoke [4].

Fire

Now that we have stable fluid simulation, we are ready to simulate fire. As discussed earlier, scalar quantities can be advected by the fluid's velocity field. We will add fuel, temperature, and smoke scalar values to our simulation. The combination of these values will be used to simulate the behavior of the fire and render it.

To keep things simple, we will assume that fuel doesn't need to be ignited. If present, fuel will give off temperature proportional to the amount of fuel present. Let's define a scalar field $\rho$ where $0 \leq \rho \leq 1$ to represent the proportion of fuel at a given coordinate in our simulation grid. The scalar field $T$ where $0 \leq T \leq T_{max}$ will represent the temperature at a given coordinate in our simulation grid. At each timestep, the new temperature $T'$ will be:

Where $T_{burn}$ is the temperature at which the fuel burns. Temperature also dissipates naturally, which can be modeled via Stefan-Boltzman Law:

The parameter $\sigma_{cool}$ is the rate of cooling. If desired, this can be tuned to be physically accurate, but this simulation is already not physically accurate. Instead, I tuned the cooling rate to get the desired look through trial and error. Also note the implicit numerical dissipation from Semi-Lagrangian advection, which will combine with the dissipation from Stefan-Boltzman law.

Other approaches opt to have a burn rate for the fuel. While this is more natural, the implicit numerical dissipation behaved well enough for my liking.

The visually simplest way to inject fuel into the simulation was as a fireball. Initially, fuel was injected as a 3D Gaussian distribution:

where $\mathbf{x}$ is the grid cell position, $\mathbf{c}$ is the center of the fireball, $r$ is the radius of the fireball, and $\sigma$ is the standard deviation controlling the spread. The raw Gaussian is then scaled by an injection rate and modulated by 3-octave FBM noise to break up the perfect symmetry and give the fireball a more organic shape.

With temperature added to our simulation, it now needs to affect how the fluid moves. The hot fluid should rise, and the cooler fluid should fall. This is our buoyancy force:

Where $\beta$ is a parameter to tune the level of buoyancy and $\textbf{j}$ is the upward vector. This term is added to the external forces term of the Navier-Stokes equation, but in practice, it will also decay to prevent the simulation from fogging up completely.

Smoke

For smoke, let's define another scalar field $S$ where $0 \leq S \leq 1$. Smoke will also be generated from fuel.

where $\gamma$ is the smoke decay rate, $k_s$ is the smoke generation rate, $\rho_{th}$ is the minimum fuel level required to produce smoke, and $\mathcal{S}$ is the smoothstep function:

The smoothstep maps fuel density smoothly from zero at $\rho = \rho_{th}$ to one at $\rho = 1$, so only fuel above the threshold contributes smoke. The threshold $\rho_{th}$ was necessary to prevent smoke from being generated at the outer edge of the fireball where fuel concentration is low, which would cause dark smoke to visually overpower the bright core of the fire.

With these scalar values in place, we can now render the fire.

Rendering and GPU Utilization

Rendering fire

At room temperature, objects emit blackbody radiation in the infrared. As temperature rises above ~700K, emission shifts into the visible spectrum, causing objects to glow — first red, then orange/yellow, then white-hot. Planck's Law describes the spectral density of light radiated by a black body at a given temperature $T$:

where $\lambda$ is the wavelength, $T$ is the temperature in Kelvin, $h$ is Planck's constant, $c$ is the speed of light, and $k_B$ is the Boltzmann constant.

So at any point in our simulation, we can sample the temperature and get the approximate spectrum of radiation that point would emit. This spectrum is responsible for making our fire simulation glow.

We can use the CIE color matching curves to translate our spectrum to sRGB values.

Luckily, the work of Wyman, Sloan, and Shirley in the Simple Analytic Approximations to the CIE XYZ Color Matching Functions allows us to translate efficiently [5].

Now that we have the ability to sample the color that the fire emits at any point in our simulation, we are now fully ready to start displaying the fire to the screen. A naive approach would be to take whatever color value is closest to the viewing perspective and render that on screen. This would show the body of the fire and how it moves, but it would lack the red to orange to yellow gradient iconic to fire because it would only be rendering the surface of the fire.

Instead, we need to turn to volumetric rendering. Specifically, ray marching will be utilized with Beer-Lambert's law. Rays are shot out from the camera's perspective, and they are iteratively marched a small distance. At each iteration, the color of the simulation is sampled, which contributes a proportion of the ray's color. These samples are added together to get the final color. This allows us to see the hot yellow inside of the fire because every bit of the fire along the ray's path contributes to the color.

GPU Utilization

As discussed earlier, the need for speed is vital to get a 3D fluid simulation to work in real time, and utilizing the GPU is essential for this requirement. The GPU was an unfamiliar domain for my coding experience before this project. I had a solid understanding of its purpose, but it was always abstracted away from me in previous use cases. I never got to interface with it with shader programs before I did this project. I used a more traditional approach to using a GPU that I would like to discuss, which had few benefits compared to the drawbacks. But, it was still a fun topic to dig into. It has led me to understand that using the GPU effectively is like taming a horse — the wrangling of race conditions, shader debugging, and rethinking algorithms for parallelism is the fuss, but once harnessed, the power is undeniable.

The Data Structure

GPUs (Graphics Processing Unit) as their name suggests were created for graphic rendering. The problem was fundamentally parallelizable - every point in a scene can be treated independently of each other. This led their designs to be as parallel as possible, but as they were developed, they became coupled with graphics. The core pipeline was to pass in the vertices of your objects to render, break those down into individual pieces called fragments in the vertex shader, and then use the fragment shader to assign color to those fragments. A common use case is to apply a texture to these fragments, so textures became a common data structure in shader programs. From a fluid simulation perspective, these concepts do not map well, but we want to harness the parallel compute.

Luckily, textures can be thought of as a more traditional and familiar data structure - an array. Fundamentally, a texture is an array of RGBA values. They hold meaning when they are rendered to the screen when they display an image, but they don't have to hold an image. Earlier, I discussed how the simulation is a grid based, and a grid is just an array with dimensions. So, textures are used to represent the simulation, and the RGBA values hold properties of the fields. The 3D vector field has its x, y, and z dimensions in the R, G, and B values respectively in the texture. Another texture holds the scalar values (fuel, temperature and smoke). With our interface into the GPU, the shader programs can implement the equations of the simulation to manipulate how the vector field and scalars change at each timestep.

The Struggles of Shader Programs

Using textures as our data structure has its pros and cons. Again, textures were originally made to apply images to objects in a scene, so writing to them wasn't a main part of their design. It is still possible with some techniques. One main challenge is that WGSL doesn't allow the same texture to be read from and written to in the same compute pass. To get around this, two textures are used to represent the same values. On every compute pass, one texture is read from, and the other texture is written to. They then alternate for each time stamp of the simulation. This process is called ping-ponging, and we use a struct to track that state.

pub struct PingPong {
  texture_a: Texture,
  texture_b: Texture,
  a_to_b: bool,
}

Alternatively, I could have explored compute shaders more, which are specifically designed for utilizing the GPU outside of graphics. These would have allowed me to use plain arrays and write less shader programs. But following GPU Gems [4], I decided to follow a traditional approach for the sake of curiosity.

Motivation

I had this project idea in the back of my mind for about a year. When I was taking a computer graphics class at UC Berkeley, my final project group had to pick a topic, one of the topics being a fire simulation based off of Andrew Chan's fire simulation project. My group took a different direction, but I've always wanted to pursue this myself. I admire Chan's way of presenting his work, which led to me to model my website off of his, build a project inspired from one of his own, and write a blog post about it following his style.

I decided to use rust here simply because I wanted to develop my rust skills, but it came equipped with the wgpu crate, which had plenty of documentation to get started. Another one of my goals was to utilize the GPU and get some experience writing programs for it.