Simulating Smoke
Antoine McNamara and Adrien Treuille
University of Washington
CSE 557 - Final Project
3/21/2002
Introduction:
For our final project, we developed a simulator which models,
animates and renders smoke according to the laws of fluid dynamics.
Animating smoke has long been a topic of interest in computer graphics,
but only in the past 5 years have researchers used physical simulations
to generate these animations. Smoke, like other fluids, is extremely
difficult to animate by hand, and there is a clear desire for simulations
that produce realistic results, preferably that allow a certain degree
of control and that are not prohibitively expensive, computationally.
We can consider slow-moving gases to be incompressible fluids, which
are modeled by the Navier-Stokes equations. Foster and Metaxes first used
the Navier-Stokes to animate gases and produced nice results on relatively
coarse grids [2]. In [4], Stam proposed a method using similar fields, but
with an implicit, semi-Lagrangian technique for updating the grids which
resulted in much simpler computations that were guaranteed to be stable.
In 2001, Fedkiw, Stam, and Jensen extended this algorithm specifically
for smoke simulation, noting that the viscosity term in the Navier-Stokes
equations didn't contribute dramatically to the motion of gases, and
proposing a technique known as "vorticity confinement" to inject velocities
back into the system that might have been lost in the numerical dissipation
of the implicit semi-Lagrangian integration scheme.
Implementation:
Our smoke simulator is divided into two parts: a physical simulator
and a renderer. During each iteration we evolve the physical simulation
by a small time step, then pass that information to the renderer. We have
created simulations and renderings both in 2D and in 3D.
i. Physical Simulation:
The central component of the physical simulator is a grid of
velocities u, representing currents of air. As in [1], we evolve
the velocity grid over time according to the Euler equation for fluid
dynamics:
Here rho is the density, p is the pressure and
f is an arbitrary force. (These equations are a simplification
of the Navier-Stokes equations that ignore viscosity.) We also wish
that the air currents preserve the mass of any smoke that passes through
them. Technically speaking, this means that the air velocity field be "divergence
free." By defining a linear operator P that projects the
velocity field onto a divergence-free velocity field, we can derive the
following equation:
The details of this derivation can be found in [4]. This equation
forms the basis for our simulation. The algorithm is as follows:
for each iteration
add the force
// f
advect the velocity along itself
// -(u del)u
project the velocity field onto a divergence-free
field // P(...)
In order to advect the velocity we use a second-order adaptive
Runge-Kutta path tracer. To compute the projection step we must solve
a Poisson differential equation. (This is a consequence of the mathematics.
We refer the reader to [4] for details.) We solve the Poisson equation
by spacially discretizing it, yielding a sparse linear system of equations.
We then use the Conjugate-Gradient method to approximate a solution.
The above algorithm gives us a field u of air currents
that evolves over time according to the laws of fluid dynamics. We then
inject two "substances" into the air: smoke and temperature. We model
these substances as a grid of scalar values. In the case of smoke, each
scalar value represents the density of smoke at a particular grid point.
For each iteration of our algorithm we advect (transport) these scalar values
through the grid along the air velocity field.
Finally, there is a feedback mechanism: as the temperature
and smoke densities get transported through the air, they affect the air
currents. Hot air rises, while dense smoke falls.
ii. Rendering:
To render the 2D animations, we simply draw the density grid
on the image plane, using OpenGL's GL_QUADs for rendering at near interactive
rates (for grids up to about 200x200). The produces very nice results,
but, of course, lack shading and self-shadowing. In 3D, we use much
more sophisticated rendering techniques.
We implemented a stripped-down ray tracer (with no reflection
or refraction rays), which casts rays through the grid and attenuates
the intensity of the ray by the amount of smoke it passes through. As
a preprocessing step, we shoot a ray from the center of each voxel toward
the light source to determine the color of the smoke at that particular
grid location (we only use a single light source, but in practice we could
do the same for multiple lights and sum the resulting intensities). We
determine the length of the path through each voxel and take the ratio over
the maximum path length as the "amount" of smoke the ray passes through.
At each voxel, the light intensity from the rest of the path is scaled by
this ratio, resulting in an exponential decrease as we proceed through
the smoke.
Once this preprocessing is complete, we cast rays from each pixel
in the viewing direction (we use a parallel projection in all of our
examples, but it trivially extends to perspective projections). Each ray
marches through the voxels in the same fashion as before, only now it uses
the path length through each voxel to decide between the preprocessed smoke
color and the color along the rest of the ray (which terminates in a background
color).
When the viewing rays exit the cube on which the grid is defined,
we use one last step to cast shadows. A ray is cast from the exit point
onto the ground plane, followed by another ray from this intersection
point in the light direction. If this ray passes through any smoke,
the background color should be attenuated to cast a shadow. We could
easily repeat the method used in the preprocessing step to accurately
determine the attenuation, but in the interest of reduced computation,
we simply test the ray for intersection with the box and then use the
precomputed value at the nearest voxel. This essentially projects the
grid onto the floor leading to aliasing of the shadow, but very little computation
is required for the shadow and the resulting shadow still serves the purpose
of grounding the smoke.
Discussion:
This project was really rewarding and lots of fun. We essentially
accomplished what we sought, but there were a few snags along the
way.
The aspect of the project that gave us the most trouble was projecting
the velocities down to a divergence-free field. In [4], the authors
suggest using a Fortran library POIS3D to solve the resulting Poisson
equation directly, but [1] suggests instead solving it as a sparse linear
system. We chose to follow the latter path, but had a hard time getting
a solver to work correctly. Initially, we weren't sure of the properties
of the matrix, and decided on the general biconjugate gradient descent
approach described in [3]. Unfortunately, we couldn't get the
linbcg solver to work under ANY conditions, even with simple 3x3
matrices. We then studied the problem more and realized our linear system
was both symmetric and positive definite. This enabled us to use a simpler
conjugate gradient solver that worked fine and produced accurate projections.
We toyed with the idea of allowing arbitrary boundaries within
the grid and having the smoke be strictly maintained within a box or
passing around objects, but never got it working. The papers discuss
these cases and suggest setting the velocity in the normal direction
of the surface to zero, but don't mention how this affects the projection
step. Clearly it must: if it didn't, setting the velocities to zero before
projecting would let the non-zero velocities bleed through the boundaries,
but setting them to zero after would mean the field was no longer divergence-free.
Therefore, the projection step must somehow be accomplished with the added
constraint that the boundaries remain zero. This might not be difficult,
but we weren't able to get it to work in the limited time we had for the
project.
It goes without saying that we encountered numerous (sometimes comical) bugs
along the way. One of the toughest to fix was that huge velocities would sometimes
sweep into the air current grid, wildly tossing the smoke about. The
final video of our "results" section illustrates this bug. We believe that
the problem was along the boundary of the grid. We solved it by setting all
velocities to zero at the boundary voxels both before and after the projection
step. This prevented any paths from being traced out of the grid during advection,
essentially creating a 1-voxel buffer between the inside and the outside
of the grid. All the simulations in the "results" section, except for the
last, use this safety buffer.
Lastly, we spent a lot of time trying to make the smoke's motion
realistic. The problem is that our simulation often looked "damp." That
is, the smoke did not twist and curl enough. While the physics do accurately
model the behavior of smoke, real smoke is never seen in a vacuum and
is always interacting with random wind and noise. Moreover, the implicit
solution methods we employed, while numerically stable, suffer from numerical
dissipation, making the problem even worse.
We addressed this issue in two ways. First, at each time step we added small,
random wind velocities. The first and second videos in our "results" section
show this. The first simulation was conducted without any random velocities,
while in the second, we initialized grid with random velocities and subsequently
injected small, random velocities at each time step. The result is much more
natural looking smoke: one does not usually experience
smoke in a perfectly still room. Our second solution was to add a vorticity
confinement force at each time step, as in [1]. This force is computed
from the velocity grid using the "curl" operator. It is consistent with the
physics in the sense that it disappears as the spacing between grid points
goes towards zero. The third and fourth videos illustrate the difference
between simulating without and with vorticity confinement, respectively.
We did not simulate vorticity confinement in 2D, as the "curl" operator
is only defined in 3D.
Results:
The following are some movies that we created during this project.
To play them you will need a
Quicktime player.
The first movie is a 2D smoke simulation. Here, we continuously inject
smoke at three points in the bottom of the frame. We do not use vorticity
confinement. The air field initially has zero velocity at all points, and
we do not add random velocity during the simulation. Therefore the smoke
rises perfectly smoothly and symmetrically.
The following video is identical to that above, except that air field
initially starts with random velocities, and we inject small random velocities
into the air field at every time step. This yields the more chaotic and natural-looking
smoke.
This video is analogous to the first video of the results section, except
in 3D. The air field initially starts out with zero velocity, and we do
not inject random velocities during the simulation. Instead of injecting
smoke at three points, we continuously inject a ball of smoke near the floor.
As in the first video, the rise of smoke is both smooth and symmetric.
This next video is like the last except that the air field starts off
with initial random velocities. Also, we have added vorticity confinement,
which makes the smoke twist and curl more naturally.
Here, we inject smoke at the bottom, slightly off-center. We also give
the smoke an initial horizontal velocity, giving the effect of a "smoke
gun." This simulation uses initial random velocities and vorticity confinement.
Finally, we present a "bloopers" video. At this point we were trying
to do object intersection with an invisible sphere. Somehow, strong velocities
enter the simulation half-way through, and the result is... well you can
see below. We go into further detail about this bug in the "discussion"
section.
References:
[1] R. Fedkiw, J. Stam, and H. W. Jensen. Visual Simulation of Smoke.
In SIGGRAPH 2001 Conference Proceedings, Annual Conference Series
, pages 15-22, August August 2001
[2] N. Foster and D. Metaxes. Modeling the Motion of a Hot, Turbulent
Gas. In SIGGRAPH 97 Conference Proceedings, Annual Conference Series
, pages 181-188, August 1997
[3] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery. Numerical
Recipes in C. Cambridge, England: Cambridge University Press,
1992.
[4] J. Stam. Stable Fluids. In SIGGRAPH 99 Conference Proceedings,
Annual Conference Series , pages 121-128, August 1999
