\documentclass[fleqn,10pt,twoside]{article}
\usepackage{parco2003}
\usepackage{epsfig}
\usepackage{graphics}

\newcommand{\ttbs}{\char'134}
\newcommand{\AmS}{{\protect\the\textfont2
  A\kern-.1667em\lower.5ex\hbox{M}\kern-.125emS}}

\title{Real-time Visualization of Smoke through Parallelizations}

\author{Torbj\o rn Vik\address[NTNU]{Norwegian University of Science and Technology,
        Department of Computer and Information Science, NO-7491 Trondheim, Norway},
        Anne C. Elster\addressmark[NTNU] and Torbj\o rn
        Hallgren\addressmark[NTNU]}

\begin{document}

% typeset front matter
\maketitle

\begin{abstract}

Animation/visualization of natural phenmena such as clouds, water,
fire and smoke has long been a topic of interest in computer
graphics, and following the increase in common processing power researchers have
started using physical simulations to generate these animations. Smoke, and other
fluids, is extremely hard to animate by hand, and simulation techniques producing
realistic results are highly applicable and desirable in areas ranging from computer
games to weather forecasts. To maintain a reasonable level of usability, these
techniques should allow a certain degree of control, and not be prohibitively
expensive computationally. However, since these tasks are extremely 
compute-intensive, powerful parallel hardware systems is needed to solve these 
tasks in real-time.

Our work combines mathematical, graphical and technical techniques 
including several parallel ization techniques in order to produce a 
{\it real-time} code that visualizes rolling and
twirling smoke in three dimensions. 
The resulting 3D parallel implementation
exploits dual cpu hardware in order to produce more than 20 frames 
per second on today's modern desktop computers.


\end{abstract}
{\bf Key words:}{\it real-time, parallel, visualizatio.n smoke, CFD, 
CG, pipelining, Intel SSE instructions}

\section{INTRODUCTION}

\label{Introduction} Early researchers in the field of visual
fluid simulations produced realistic results by applying
mathematically generated textures to simple geometric structures
\cite{Gard85,Perl85}. This approach did not count for interaction
with solid objects and left the animators with the tedious work of
parameter tweaking.

The models and algorithms presented in
\cite{Fedk01,Stam99,Fost97}, all demonstrate a clear evolution
within visual simulations of fluids: In 1997, Foster \& Metaxas
published their paper, \textit{Modeling the Motion of a Hot,
Turbulent Gas} \cite{Fost97}, introducing an CFD algorithm for
visual fluid simulations. They extracted the parts of
Navier-Stokes' equations that affect the interesting visual
features of smoke, such as convection and vorticity, thus creating
a model capable of simulating realistic smoke on relatively coarse
grids. Their model was only stable for small time steps because of
an explicit integration scheme; large velocities required small
timesteps, and the result was animations evolving too slowly. Stam
was the first to present a stable algorithm solving the full
Navier-Stokes equations \cite{Stam99}, independent of the time
steps and implemented a version that "allows an animator to design
fluid-like motions in real time". This was achieved using a
semi-Lagrangian approach together with implicit solvers, allowing
for much larger time steps and a faster evolving simulation.
Fedkiw, Stam and Jensen \cite{Fedk01} pointed out that the
first-order integration scheme used by Stam caused too much
numerical dissipation and small-scale details just tended to die
out. Fedkiw et al. presented a model aiming for visual simulation
of smoke, compared to the general fluid simulator of Stam. Their
model was based on the work of Foster \& Metaxas \cite{Fost97} and
Stam \cite{Stam99}, but featured some modifications for improved
visual quality. They reported frame rates ranging from 10 fps to 1
fps for their fastest, low-quality simulations.

For a real-time environment aiming for a visual imitation of the
real world, these frame rates are, however, still not sufficient. We present
a configuration of parameters, sub-algorithms and numerical
solvers, all analyzed with respect to visual quality,
computational efficiency and parallelism, eventually aiming for
frame rates of 20 and above on today's modern desktop hardware.

\section{SIMULATION MODEL}

\label{SimulationModel}Our simulation model is based on the one
presented by Fedkiw et al. \cite{Fedk01}.
The algorithm dices the computational doimain into uniform
sub-cubes known as \it{voxels} and uses two of copies these 3D grids
-- one grid contain ing the data for the previous frame while 
the other stores the data enerated for the current frame.
Each grid defines the velocity vector \textit{u},
and values such as density \begin{math}\rho\end{math}, temperature
\textit{T}, and user-defined forces (e.g. wind).\\

Our smoke simulation algorithm consists of 5 numerical steps
followed by a step for calculating the lighting and a step 
for rendering giving us the
following 7-step algorithm for each frame:

\begin{enumerate}
\item{\bf{Add Forces}, including the bouyancy force (gravity and effects of hot air rising)
 and user-defined forces}, and update the velocites by scaling the forces with
the timespet.
\item{\bf{Self-Advect velocities} through Euler's equations 
(assuming  inviscid fluid)}:
\begin{equation}
\frac{u-u*}{\Delta t} = -\left ( u \cdot \nabla \right )u + f
\label{eqSolveNS2}
\end{equation}

\item{\bf{Solve Poisson's Equation} for the pressure using CG
\begin{enumerate}
	\item{Build Gradient vectors}
	\item{CG}
\end{enumerate}
\item{\bf{Update (project) temporary velocities} by subtracting
the pressure gradient scaled by the time step in order
to conserve mass (assuming incompressible fluid).}

\item{\bf{Advect scalar values}, i.e.  update density $\rho$and temperature $T$
with respect to the domain similarly to step 2}
\item{Calculate lighting}
\item{Render image}
\end{enumerate}

The calculation of the next frame is dependent of the results from
Step 6. As a result, only Step 7 -- the rendering step -- may
overlap the computations with the next frame. This implies that
the simulation algorithm itself is inherently serial. However,
each computational step, except rendering (which may overlap with
the calculations of the next frame), may benefit from
parallelization as will be shown in the next section.


\section{SERIAL OPTIMIZATIONS}
In order to get real-time performance we first optimized
our serial code in order to have a good starting point
for our code.

We first built a 2D prototype that we analyzed and profiled 
in order to pinpoint possible processing hotspots in
the algorithm. See Figure ****  Here, velocities were defined at voxel
centers as in Fedkiw et al. \cite{Fedk01}. Notice that the two advection
steps occupied about 90\% of the CPU time. 

Our 3D implementation adopted the approach of Stam\cite{Stam99} and
Foster et al.\cite{Fost97}, with velocities defined at voxel
centers, thus reducing the number of linearly interpolated
off-grid values. This also helped our parallel implementation.

\subsection{Possoin's Solver}
Solving Poisson's equation requires a linear solver handling
sparse Symmetric Positive Definite (SPD) matrices in an efficient
manner, e.g. relaxation methods \cite{Fost97,TracXX} or iterative
solvers \cite{Kincaid**, Young******Fedk01,Stam99,Henk94,Shew94}. 

Foster et. al uses a incomplete Cholesky preconditioner as part 
of their Conjugate Gradients (CG) solver \cite{Shew}. 
However, since we know that for real-time
smoke simulations, the time-steps and changes between each
frame are realtively small, we skip the pre-conditioner and
use the solution from the previous frame as an initial guess.
As a result, our code only uses 5 iterations to converge.
 


\section{PARALLELIZATION TECHNIQUES}\label{Parallelization}

The computational  steps described in the previous section,
consist of a series of iterations over all the voxels
in the domain. Iterative operations are usually very suited for
parallelism. Step 3 (Solve Poisson's equation) is divided into 
two individual steps because of the
synchronization point in between the two iterative constructs
which is of importance when parallizing the algorithm.

Build-gradient vector and CG. This step can still be run in
parallel locally, by dividing the iterative constructs and
assigning one part to each node.

The rendering is performed by constructing a set of transparent
planes from the volume data, using OpenGL [ 12] as software API to
graphics hardware.


We have implemented one serial and two parallel versions
of the simulation model: one data-parallel and one modified data
parallel version, where the latter supports the technique of
pipelining.

\subsection{Data Parallel Implementation}\label{DataParallelImplementation}

The data-parallel version divides the different iterative
constructs of the algorithm into smaller parts and distributes
them to the different nodes, and performs administration and
synchronization using Win32 API calls. This of course adds some
overhead, but when running the simulation on relatively small
number of CPUs (2 - 4), this overhead should be negligible.

Section \ref{ResultsDataParallel} presents benchmarks and analysis
of the data-parallel approach. In short, it proves a significant
speedup by exploiting parallelism but also uncovers that one of
the CPUs is left idle while the other is running the OpenGL
rendering code.

\subsection{Pipelined Implementation}\label{PipelinedImplementation}

The dataflow and dependencies within a single frame of animation
show that a frame does not require the entire previous frame to
finish before it can start processing. The data-parallel
implementation also uncovered that certain serial parts of the
algorithm leaves one or more CPUs idle, and, together with the
partial frame-dependencies, this opens for the technique of
pipelining: jobs from two different frames can be interleaved, as
their dependencies do not require them to be processed in exact
sequential order.

We have implemented a pipelined version exploiting this situation,
using a round-robin queue to distribute jobs among the worker
threads and administrate the dependencies. Each iterative
construct of the algorithm is divided into smaller parts, similar
to the data parallel version. Each job is set up with a set of
pointers to the jobs it opens for. Each job is also set up with an
integer telling the number of parent jobs it has, and how many of
the parents that are currently actually finished. Whenever a job
is completed, each of its children gets their counter increased.
Whenever a job's counter equals its number of parents, it is
available for processing. A job's counter is reset immediately
after a thread has started working at it. This way, while one
thread performs the rendering, other threads may start processing
the next frame and complete it either partially or completely.

The efficiency of this technique depends on the performance of the
graphics card, or to be more specific, how much time that is spent
performing the rendering compared to the time spent performing
numerical simulation. The efficiency increases as the time
required for rendering and numerical processing becomes more and
more the same. For the pipelined approach to produce any speedup
at all, there must be a reasonable amount of processing time that
is unconditionally serial. If the rendering step completes too
fast, the other CPU will not be able to process anything of the
following frame. In this case, the pipelined approach will
probably slow down the system due to the overhead added through
the queueing system.

\subsection{Vectorization on Intel CPUs}
Vectorization of iterative program constructs is a technique
usually employed in the context of the so-called vector computers.
However, the current Intel family of processors includes a set of
SIMD instructions, known as SSE (Streaming SIMD Extension) which
work on floating point values in groups of four, theoretically
performing four times as many calculations per clock cycle on a
single processor by utilizing the multiple arithmetic unit
available on today's Intel CPUs.

We modified most of the simulator to exploit these SSE
instructions, and experienced about 10\% average speedup. Some
parts of the algorithm are more suited for vectorization than
others, especially the method of Conjugate Gradients. The more
complex steps, such as the advections, require interpolation of
off-grid values using eight different on-grid values not placed
contiguously in memory. An SSE version of the advection steps
hence require the values to be gathered into special vector data
types in groups of four. This gathering and conversion adds
overhead to the code, thus reducing the speedup.

SIMD instructions can increase the memory bus usage significantly
compared to standard floating-point operations since the bus must
deliver more data to keep up with the parallel computations. In
other words, the memory bus may end up being an even tighter
bottleneck than before.

\subsection{Parallelization of lighting calculations}
The lighting calculations cannot be treated as ordinary iterative
operations, since the order of traversal (of the voxels) depends
on the light's direction and where the rays enter the volume. We
parallelize Fedkiw et al's lighting model \cite{Fedk01} by
grouping the light rays into distinct sets, and assigning each set
to a thread. The model supports parallel light rays entering the
volume from an arbitrary angle.

\section{RESULTS}\label{RESULTS}

\subsection{Results for parallel implementations}\label{ResultsDataParallel}

Figure \ref{figure6} shows the speedup gained by our data-parallel
implementation, compared to the serial version and pipelined
version.
\begin{figure}
\begin{center}
\psfig{file=parco_figure5.eps} \caption{The speedup of data
parallel and pipelined implementation vs. serial.}\label{figure6}
\end{center}
\end{figure}
The chart shows that there is indeed a speedup with the data
parallel implementation, varying from 20\% to 60\% depending on
the simulation parameters. Still, 60\% speedup when moving from
one to two CPUs is far from linear scaling. A partial explanation
is that one thread is performing rendering while the other threads
stall for some time, waiting for the first to finish. This is
confirmed by the numbers reported by the MS Task Manager, giving a
CPU utilization ranging from 60\% for coarse grids to 90\% for
fine grids. However, even with rendering disabled, the speedup is
still not higher than about 66\%, indicating that there are other,
more substantial, explanations. Some of the processing power is of
course lost due to overhead added when our code administrates and
distributes tasks among the worker threads. But, as stated in
section \ref{DataParallelImplementation}, this overhead should be
close to zero compared to the different calculations involved in
the algorithm. Most likely, the poor scaling is caused by race
conditions for a shared resource somewhere in the system, and
tests point out the memory bus as being the probable cause. The
data required by a CPU to fulfill its work is way larger than its
local cache, and the memory bus is unable to deliver data fast
enough for the program to scale linearly.

Figure \ref{figure6} also illustrates the evident speedup gained
by the pipelined approach, compared to the data-parallel, but the
efficiency is most likely still limited by memory bus bandwidth.

\subsection{Visuals}\label{Visuals}

\begin{figure}
\begin{center}
\psfig{file=parco_figure6.eps} \caption{Frames [60, 110, 160] of a
smoke simulation, domain size 16x16x64, 5 CG iterations, 40 fps on
[DualFAST].}\label{figure7}
\end{center}
\end{figure}\begin{figure}
\begin{center}
\psfig{file=parco_figure7.eps} \caption{Frames [60, 110, 160] of a
smoke simulation, grid sized 32x32x128, 5 CG iterations, fps =
4.1.}\label{figure8}
\end{center}
\end{figure}\begin{figure}
\begin{center}
\psfig{file=parco_figure10.eps} \caption{Frame 60, 110 and 160 of
a smoke simulation, grid sized 32x32x128, 5 CG iterations,  fps =
4.1. Random forces.}\label{figure10}
\end{center}
\end{figure}
Figure \ref{figure7} presents a series of pictures from our
implementation. It runs about 40 frames per second on a Dual Intel
Xeon 1.7 GHz, simulating 16 384 voxels. The smoke moves in a
realistic fashion and interacts correctly with the obstacle
present in the domain, but suffers from the numerical dissipation
and lack of detail pointed out by Fedkiw et al. \cite{Fedk01}.

Figure \ref{figure8} presents a simulation with the exact same
parameters as Figure \ref{figure7}, but features eight times as
many voxels (131 072). The detail is greatly improved, but the
result is only 4 frames per second.

Figure \ref{figure10} illustrates the effect of introducing small
random forces into the domain, and show, in comparison with Figure
\ref{figure8}, that the smoke looks less "smooth" and for some
situations more realistic.

\section{CONCLUSION}\label{CONCLUSION}

The results showed that it is possible to run an optimized
implementation of the algorithm outlined by Fedkiw et
al.\cite{Fedk01} that yields 5 to 40 frames per second on today's
modern home computing hardware. Algorithmic optimizations,
parameter tuning and efficient programming techniques made the
implementation capable of animating smoke, still with a
satisfactory level of realism and speed in real-time.

We have involved two types of optimizations: model and code. The
former modified the algorithm and sacrificed accuracy for speed,
while the latter simply increased the efficiency of the program.
Three important modifications were done on the model from [1]:
\begin{samepage}
\begin{itemize}\item[-]{\textit{Defining the
velocity at voxel centers instead of voxel
faces}}\item[-]{\textit{Fast matrix-free CG (no
preconditioner)}}\item[-]{\textit{Leaving out the vorticity
confinement}} \end{itemize}\end{samepage}

These three decisions unquestionably reduced the visual quality of
the smoke, a fact supported by a comparison of the visual results
presented by Fedkiw et al.\cite{Fedk01} with the visual results of
this our implementation (see section \ref{Visuals}). Nevertheless,
the modifications also reduced the computational complexity of the
algorithm, thus making it possible to calculate significantly more
frames per second. It is interesting to notice that the number of
CG iterations made little or no difference on the visual
simulations; a surprising result when seen in the light of earlier
discussions \cite{Fedk01}. This effect should be investigated
further.

The speedups gained by parallel implementations of the algorithm
proved it very suitable for parallelism, but also revealed the
memory bus as being a limitation for memory intensive tasks. The
pipelined approach presented in \ref{PipelinedImplementation}
introduced a technique for dynamic job distribution and
pipelining, thus reducing the penalty caused by serial sub-tasks
such as rendering. The technique of pipelining is well-known,
especially within computer hardware engineering, but our results
showed its relevance within the field of CFD visualization. And as
most of today's SMP systems still share a common graphics card,
the technique could apply to a wide range of visualization
software. Results also showed that certain parts of the algorithm
gained significant speedups as a result of Intel's SSE
instructions, especially the method of Conjugate Gradients.

Whether it is possible or not to run a visual CFD simulation in
real-time, depends on the desired visual accuracy and quality.
Higher software and algorithmic efficiency will allow for higher
accuracy on a on a given hardware platform, but the numerical
simulations will always have a potential for improvement.

\begin{thebibliography}{9}
\bibitem{Fedk01} R. Fedkiw, J. Stam, and H. W. Jensen. Visual Simulation of Smoke.
                 In SIGGRAPH 2001 Conference Proceedings, Annual Conference Series , pages 15-22, August 2001
\bibitem{Stam99} J. Stam. Stable Fluids.
                 In SIGGRAPH 99 Conference Proceedings, Annual Conference Series, pages 121-128, August 1999.
\bibitem{Fost97} N. Foster and D. Metaxas. Modeling the Motion of a Hot, Turbulent Gas.
                 In SIGGRAPH 97 Conference Proceedings, Annual Conference Series, pages 181-188, August 1997.
\bibitem{Henk94} Henk A. Van Der Vorst. Lecture Notes on iterative methods.
                 Mathematical Institute, University of Utrecht, Budapestlaan, the Netherlands. June 4, 1994.
\bibitem{Shew94} An Introduction to the Conjugate Gradient Method Without the Agonizing Pain Edition 1.25.
                 Jonathan Richard Shewchuk, August 4., 1994, http://www-2.cs.cmu.edu/~jrs/jrspapers.html
\bibitem{BartXX} Peter Bartello and Stephen J. Thomas. The Cost-Effectiveness of Semi-Lagrangian Advection.
                 Dr. P. Bartello, RPN, 2121, voie de Service nord, Route Transcanadienne, Dorval (Quebec) H9P 1J3.
\bibitem{TracXX} Fred T. Tracy. A Comparison of a Relaxation Solver and a Preconditioned Conjugate Gradient Solver in Parallel Finite Element Groundwater Computations.
                 Engineer Research and Development Center, Major Shared Resource Center. http://www.hpcmo.hpc.mil/Htdocs/UGC/UGC01/paper/fred\_tracy\_paper.pdf
\bibitem{Gard85} G. Y. Gardner. Visual Simulation of Clouds.
                 Computer Graphics (SIGGRAPH 85 Conference Proceedings), 19(3):297-384, July 1985.
\bibitem{Stei98} C. M. Stein, N L Max. A Particle-Based Model for Water Simulation.
                 (SIGGRAPH 98 Conference, july 19-24).
\bibitem{Perl85} K. Perlin. An Image Synthesizer.
                 Computer Graphics (SIGGRAPH85 Conference Proceedings), 19(3):287-296, July 1985.
\bibitem{OpenMP} OpenMP,
                 http://www.openmp.org
\bibitem{OpenGL} OpenGL,
                 http://www.opengl.org
\bibitem{Vik03}  Real-time visual simulation of Smoke, Master thesis at
NTNU,\\
                 http://www.tihlde.org/\~\ torbjorv/diplom15.doc
\end{thebibliography}

The simulation model used is based on the assumption that smoke
can be modeled as an inviscid, incompressible fluid, using Euler's
equations. The algorithms are largely based on the ones presented
by Fedkiw, Stam and Jensen in their paper entitled \textit{Visual
Simulation of Smoke}, and are analyzed with respect to efficiency
and visual quality. All the algorithms involved are discussed from
a parallel programming point of view, and two different parallel
SMP implementations are presented. More specifically, Intel based
SMP systems running MS Windows were used as test-beds.

We introduce techniques for dynamic job distribution and
pipelining, thus reducing the penalty caused by serial sub-tasks
such as rendering. The technique of pipelining is well-known,
especially within computer hardware engineering, but this work
shows its relevance within the field of visualization. And as most
of today's SMP systems still share a common graphics card, the
techniques could apply to a wide range of visualization software.
Results showing that certain parts of the algorithm gained
significant speedups by using Intel's SSE instructions, especially
the Conjugate Gradients method, are also presented.

The simulation algorithm consists of five main numerical steps:\\
\begin{enumerate}
\item{
\textit{\textbf{Add forces}}\\
Initially, the velocity field is updated with respect to all
forces. Both the buoyancy force (\ref{eqBuoy}) and any
user-defined forces are multiplied by the time step
\begin{math}\Delta t\end{math}, added up with the velocity from
the previous frame, and stored in the grid for the current frame
(\ref{eqAddForces}).
\begin{equation}
u = u + \Delta t*\left (f_{buoy} + f_{user defined}\right )
\label{eqAddForces}, \roman{where}
f_{buoy} = -\alpha\rho + \beta(T-T_{amb}) \label{eqBuoy}
\end{equation}
The buoyancy force works only in the vertical direction and
accounts for both gravity,\begin{math}-\rho\end{math}, and the
effect of hot air rising, \begin{math}T - T_{amb}\end{math}. Two
constants, \begin{math}\alpha\end{math} and
\begin{math}\beta\end{math}, should be set up with physically realistic
values.}

\item{ \textit{\textbf{Self-Advect Velocities}}\\ The smoke's
velocity is described by Euler's equations, assuming that smoke
can be modeled as an incompressible, inviscid fluid.
\begin{equation}
\nabla u = 0 \label{eqNS1}
\end{equation}
\begin{equation}
\frac{\delta u}{\delta t} = -\left ( u\cdot\nabla\right )u-\nabla
p + f\label{eqNS2}
\end{equation}
where f is the sum of the buoyance force (\ref{eqBuoy}) and all user
defined forces. The variable $\nabla p$ describes the gradient
of the pressure at each grid point (voxel). Since we are providing a visualization
of the phenomea rather than computing actual physics, we can subgrids of 
size 1, saving a lot of scaling factors. E.g the gradient pressure in the x-direction
at voxel [i, j, k] is defined as:
\begin{eqnarray}
\nabla p_{i, j, k}^x = p_{i+1, j, k} - p_{i-1, j, k}\nonumber\\
\label{eqPGradient}
\end{eqnarray}
Similarly, the velocity gradient in the x-direction is defined as
\begin{eqnarray}
\nabla u_{i, j, k}^u = u_{i+1, j, k} - u_{i-1, j, k}\nonumber\\
\label{eqUGradient}
\end{eqnarray}
Similar equations are used for the y and z directions.

According to Euler's equations, both mass (\ref{eqNS1}) and
momentum (\ref{eqNS2}) should be conserved. We solve equation
(\ref{eqNS2}) without the pressure term over a time step
\begin{math}\Delta t\end{math}, as \cite{Fedk01}:
\begin{equation}
\frac{u-u*}{\Delta t} = -\left ( u \cdot \nabla \right )u + f
\label{eqSolveNS2}
\end{equation}
The advection operator \begin{math}-\left ( u \cdot \nabla
\right)u\end{math} is solved using a semi-Lagrangian method
\cite{Stam99}: A departure point \textit{x*} is computed for each
voxel (i, j, k) at position \textit{x}, using the trapezoidal
method \cite{BartXX} as an implicit integration scheme:
\begin{equation}
x* = x - \frac{\Delta t}{2}\left (u\left (x,t \right ) + u \left
(x - \Delta t u \left (x,t\right ),t - \Delta t \right ) \right )
\label{eqTrapezoidal}
\end{equation}
\textit{t} represents the time for the current frame of animation,
\begin{math}\Delta t\end{math} is the time step and \begin{math}u \left (x, t\right )\end{math} gives the velocity at position
\textit{x} and time \textit{t}. The velocity is then transferred
from the departure point to a intermediate velocity \textit{u*} at
the arrival point (the current voxel):
\begin{equation}
u* = u \left (x*, t \right ) \label{eqTransferVelocity}
\end{equation}
Whenever \textit{x*} is an off-grid position, the velocity
\textit{u} is computed using linear interpolation.}

\item{\textit{\textbf{Solve Poisson's Equation}}\\
Next, Poisson's equation for the pressure (Eq 10) is solved.
\begin{equation}
\nabla ^2 p = \frac{1}{\Delta t}\nabla u * \label{eqPoisson}
\end{equation}
Solving Poisson's equation requires a linear solver handling
sparse Symmetric Positive Definite (SPD) matrices in an efficient
matter, e.g. relaxation methods \cite{Fost97,TracXX} or iterative
solvers \cite{Fedk01,Stam99,Henk94,Shew94}. Our solver uses a
matrix-free, unpreconditioned version of the method of Conjugate
Gradients \cite{Shew94} and starts with an initial guess equal to
the solution from the previous frame.}

\item{\textit{\textbf{4. Project Velocities}}\\
The intermediate velocity \begin{math}u*\end{math} from step 2 is
forced incompressible (projected) by subtracting the gradient of
the pressure
\begin{equation}
u=u*-\Delta t \nabla p \label{eqProject}
\end{equation}}
\item{\textit{\textbf{Advect Scalar Values}}\\
The smoke's temperature and density need to be moved (advected)
through the domain in a similar manner as the velocity in step
2.
\begin{eqnarray}
\frac{\delta T}{\delta t} = -\left ( u \cdot\nabla \right )T\\
\frac{\delta \rho}{\delta t} = -\left ( u \cdot\nabla \right )\rho
\end{eqnarray}
A departure point \textit{x*} is calculated for each voxel and the
values (\textit{T} and \textit{r}) are transferred to the arrival
point (the voxel's center).}
\end{enumerate}

This five-step algorithm produces a field of scalar values for
each time step, representing the smoke's density throughout the
computational domain. The algorithm has a computational complexity
scaling linearly with the grid's resolution (number of
grid-points).
\end{document}
