CGI2014 manuscript No.
(will be inserted by the editor)

An Optimization Approach for Designing Fluid Flow Fields
Syuhei Sato · Yoshinori Dobashi · Kei Iwasaki · Hiroyuki Ochiai ·Tsuyoshi Yamamoto · Tomoyuki Nishita

Abstract In entertainment applications, such as computer games and movies, animators are often requested
to synthesize a realistic fluid flow with a particular behavior. This paper presents a method to help animators
meet such requirements. Our method enables the user
to design a realistic fluid flow. In order to generate thedesired fluid flow in a realistic form, the flow is calculated by solving a minimization problem subject touser-specified constraints and the Navier-Stokes equations. The minimization problem is eﬃciently solved byrepresenting the flow field by a linear combination of orthogonal basis functions. Consequently, the fluid flow isobtained by simple matrix operations.
Keywords fluid simulation · incompressible NavierStokes equation · Laplacian eigenfunctions · design
1 Introduction
Visual simulation of fluids has become one of the mostimportant research topics in computer graphics. Many
methods have been proposed for simulating smoke, water, fire, and so on [1]. Most of the recent methods are
based on computational fluid dynamics to create realistic animations and these are used in many applications
such as movies and computer games. However, one ofthe problems is the expensive computational cost. In
Syuhei Sato, Tsuyoshi Yamamoto
Hokkaido University
Yoshinori Dobashi
Hokkaido University, JST CREST, UEI Research
Kei Iwasaki
Wakayama University, UEI Research
Hiroyuki Ochiai
Kyushu University, JST CREST
Tomoyuki Nishita
UEI Research, Hiroshima Shudo University

those entertainment applications, rendering fluid witha particular motion is often requested. Animators usually attempt to create the desired motion by repeating
fluid simulations with diﬀerent parameter settings untila satisfactory result is obtained. However, this is an extremely tedious and time-consuming task incurring an
expensive computational cost.
Many methods have been proposed to address thisproblem. Controlling the fluid simulation is a promising
approach [2, 3, 5]. However, a trial-and-error process isstill required to tune the control parameters.
Several recent researches have focused on a diﬀerent approach: using low-resolution fluid simulation as
a guide to produce high-resolution results. In this approach, low-resolution fluid simulation is used for designing the fluid motion, which is subsequently converted into a high-resolution animation. Some methods control the simulation using low-resolution simulation results [4, 7]. These methods allow the user toeﬃciently generate the desired high-resolution animation. However, creating a particular fluid motion is still
a diﬃcult and tedious task even when low-resolution
fluid simulation is used.
To address this problem, we propose a method thatenables the user to design a realistic fluid flow. Our
method is based on grid-based fluid simulation and issuitable for synthesizing animations of incompressible
flow such as smoke and fire. The key concept behind our
method is to solve a minimization problem subject tothe user’s constraints and the Navier-Stokes equations.
Our method decomposes the flow field into low andhigh frequency components. The minimization problem
is solved for the low frequency components to control
the overall motion, and the high frequency components
are subsequently added to create detailed motion. The
low frequency components of flow fields are represented

2Syuhei Sato et al.

(our method)
renderervu(y, tn)
uservl(y, tn)v(x, tn)minimizer +Fvh(x, tn)
Duser constraints

vns(x, tn)

NS solverv(x, tn-1)

Fig. 1 Overview of our system. D is a delay buﬀer thattakes a flow field as input and outputs it at the next time
step. F is a low pass filter.

calculated on the high-resolution grid whose grid size isNh . Then, vns (x, tn ) is decomposed into low and highfrequency components. vh (x, tn ) is obtained as high frequency components of vns (x, tn ). Our system computes
an incompressible flow field vl (y, tn ) by solving a minimization problem (see Section 3). Our minimizer takes
vu (y, tn ) and [F ∗ vns ](y, tn ) as input, where F is alow pass filter, and ∗ is the convolution operator. After computing v(x, tn ), scalar quantities, such as smoke
density, are advected and visualized by using a volume
rendering technique. The user then modifies the constraints to create the desired visual eﬀect.
3 Definition of minimization problem

by a linear combination of incompressible orthogonalbasis functions, which ensures that the velocity fieldsatisfies the divergence-free condition. Moreover, ourmethod uses 2D basis functions for 3D flow fields, andthe fluid flow is generated at a low computational cost.
Since the creation of the flow field is formulated as aleast squares minimization problem using orthogonalbasis functions, we can compute the flow field using
simple matrix operations.
2 Overview of our system
Fig.1 shows an overview of our system. The flow fieldv is calculated at discrete time steps tn (n = 0, 1, · · ·),where the time interval is ∆t. We decompose the flow
field v into low and high frequency components.
v(tn ) = vl (tn ) + vh (tn ),(1)where vl and vh are low and high frequency componentsof the flow field v, respectively. In our computation, v
and vh are sampled at each grid point x, where the gridsize is Nh . vl is sampled at each grid point y, wherethe grid size is Nl that is smaller than Nh . Our method
calculates desired fluid flow vl in low frequencies, andthen adds high-frequency components vh to vl . Details
are described in the following.
In order to generate the desired flow field, the user
first places a set of control points at arbitrary positions as indicated by the blue spheres in Fig.1. The user
then specifies the velocities at the control points (thered arrows). The user-specified velocity of the l-th control point at time step tn is denoted as su (yl , tn )(l =0, 1, · · · , Nu − 1), where yl is the position of the l-th
control point and Nu is the number of control points.
Then, our system interpolates su (yl , tn ) in order to obtain a spatially continuous flow field vu (y, tn ) by using
radial basis functions. On the other hand, a simulated
flow field vns (tn ) is obtained by solving the NavierStokes equations using the flow field v(tn − 1) computed at the previous time step (see Fig.1). vns (tn ) isIn our method, the low frequency components vl (y, tn )of the flow field are represented by a linear combinationof basis functions, Φi , as follows,vl (y, tn ) =N−1∑

wi (tn )Φi (y),(2)i=0where wi (tn ) is the coeﬃcient for the i-th basis function and N is the number of basis functions. We use
Laplacian eigenfunctions [6] for Φi (y). This enforcesthe divergence-free condition on the flow field. For reducing computation and storage costs, 3D flow fields
are represented by using 2D basis functions (see Section 5). Our task is then to determine the coeﬃcientw(tn ) = (w0 (tn ), w1 (tn ), · · · , wN −1 (tn )) at each time
step so that a realistic flow with the desired behavior
is generated. The user can specify the constraint velocities at multiple positions in space and time in order todesign the desired flow field. The Navier-Stokes equations are taken into account in computing the flow. The
coeﬃcient w(tn ) is obtained by minimizing the sum of
two error functions, Eusr and Ens . That is,arg min(Eusr (tn ) + αEns (tn )),(3)w(tn )

where α is a user-specified constant used to adjust theinfluence of the two error functions described below.
Eusr (tn ) measures the diﬀerence between the userspecified velocities and the resulting flow field. Eusr (tn )is defined by the following equation.
Eusr (tn ) = ||vu (y, tn ) −N−1∑

wi (tn )Φi (y)||2 ,(4)i=0where || · || is the L2 norm of the function. vu (y, tn )represents the constraint flow field and is defined byinterpolating user-specified velocities su (yl , tn ) by using
radial basis functions. That is,vu (y, tn ) =N∑u −1
l=0

bl (tn )φ(|y − yl |),(5)An Optimization Approach for Designing Fluid Flow Fields

z3

computed by the above equation is zero if i ̸= j. Therefore, A is a diagonal matrix and its diagonal elements
can be precomputed. Next, i-th element ci of c(tn ) is
given by ci = {vu (y, tn )+α[F ∗vns ](y, tn )}·Φi (y). Our
system computes the dot product between functions in
the above equations on a grid with a specified resolution Nl . Then, the coeﬃcient wi (tn ) is computed bywi (tn ) = ci /aii . The above approach allows us to compute the coeﬃcient vector w(tn ) very eﬃciently since
we do not have to compute A at run time. All we needto do at each time step is to project vu + α[F ∗ vns ]onto Φi , and compute wi (tn ) = ci /aii .

yxFig. 2 Slices for representing a 3D flow field.

5 Basis functions
where φ is a radial basis function and bl (tn ) is obtainedby solving the following linear system of equations.
N∑u −1bl (tn )φ(|ym − yl |) = su (ym , tn ),

(6)

l=0

where m = 0, 1, · · · , Nu − 1. For φ, we use Gaussianbasis functions.
Ens (tn ) measures the diﬀerence between the flow
field and the incompressible inviscid Navier-Stokes equations. By discretizing the incompressible inviscid NavierStokes equations with the time step ∆t(= tn − tn−1 ),we obtain the following equation.
u(tn ) = G(u(tn−1 )),(7)
1G(u(t)) = u(t) + ∆t{−(u(t) · ∇)u(t) − ∇p + f }. (8)ρLet us now assume that we have obtained the flow fieldv(x, tn−1 ) at time tn−1 . In order to satisfy the NavierStokes equations, the flow field at time tn needs to becomputed according to Eq.(7). That is, vns (x, tn ) =G(v(x, tn−1 )). We use low frequency components of vns
for defining Ens . We define Ens (tn ) as the diﬀerence
between vl (y, tn ) and low frequency components of vns :Ens (tn ) = ||[F ∗ vns ](y, tn ) −N−1
∑wi (tn )Φi (y)||2 .

(9)

i=0

Finally, vh (x, tn ) is obtained as follows: vh (x, tn ) =vns (x, tn ) − [F ∗ vns ](x, tn ).
4 Solving minimization problem
Given the definitions of the error functions in the previous section, we can now solve the minimization problem
given in Eq.(3). By taking a derivative of Eq.(3) withrespect to the coeﬃcient wi (tn ), we obtain the matrix
equation Aw(tn ) = c(tn ), where A is an N ×N matrix,w(tn ) and c(tn ) are N dimensional column vectors. The(i, j)-th element aij of A is given by aij = (1+α)(Φi (y)·Φj (y)), where · indicates the dot product between twofunctions. Since Φi are orthogonal basis functions, aij

As we have described previously, the Laplacian eigenfunctions are used as the basis functions. In our method,
we use basis functions defined on a 2D rectangular gridfor a 3D flow field. Although 3D basis functions exist[6], we found that those functions are not suitable forour purpose. In order to design a 3D flow field, we needa large number of basis functions. We tried 512 basis
functions but we found that they were not suﬃcient.
Furthermore, the storage cost for 3D basis functions
sampled on a 3D grid is very expensive. We can compute the basis functions on the fly but we found that
this significantly slows down the computation.
To address these problems, we propose to use the2D basis functions for the 3D flow fields. For using 2D
basis functions, we subdivide the simulation space using
three sets of slices, indicated in red, green, and blue in
Fig.2. The red, green, and blue slices are orthogonal tothe x, y, and z axes, respectively. On the red slices, theyz components of the flow field are represented by 2Dbasis functions. Similarly, for the green and blue slices,
the zx and xy components, respectively, are representedp
by 2D basis functions. Let us define the flow field vyz,rqvzx , and vxy on the p-th red, q-th green, and r-th blueslices, respectively. The 3D flow v at grid point (p, q, r)
is generated by combining these three flow fields:
1 pqrvl (p, q, r) = (vyz(q, r) + vzx(r, p) + vxy(p, q)). (10)2Note that the 3D flow field computed in this way satisfies the divergence-free condition. This method significantly reduces the storage cost for the basis functions
sampled at the grid points and the pre-computation
time for calculating A, because slices with the same
size can share the same basis functions. If the shapeof the simulation space is a cube, we can use the same
basis function for all the slices.
6 Results
This section shows some examples created using ourmethod. We used a desktop PC with a CPU (Intel4Syuhei Sato et al.

Fig. 3 Example of fire.

Core i7 2600K, 16GB memory) and GPU (NVIDIAGeForce GTX TITAN) to compute all the examples
shown in this section. In all examples, the number ofbasis functions (N ) is 256, the low resolution and highresolution grid sizes (Nl and Nh ) are 32 × 32 × 48 and256 × 256 × 384, respectively. The computation time for
updating the flow field is 0.7sec/frame. The videos corresponding to the following examples can be found inthe supplemental material.
Fig.3 shows an example of a fire animation with acircular swirling motion that is simulated by setting
five disk-shaped heat sources at the bottom of the simulation space. The left-upper inset of Fig.3 shows thecontrol points. Our method can be used to render fire
simulations with buoyancy forces.
Next, Fig.4 shows an animation of a tornado. The
right-bottom inset of Fig.4 shows the control points.
By moving the control points, the user can design themotion of the tornado.
7 Conclusion
We have proposed a method for designing incompressible and inviscid flow fields. The user can specify aset of constraints in order to create the desired flow
field. The method computes realistic flow fields by solving the minimization problem such that user-specified
constraints are satisfied. The Navier-Stokes equations
are also taken into account in solving the minimization problem. The low frequency components of flow
fields are represented by incompressible orthogonal basis functions and this ensures the divergence-free conditions of the flow. We also proposed a method for representing a three-dimensional flow field using a set ofFig. 4 Example of tornado.

two-dimensional basis functions. We demonstrated thecapabilities of our method with a set of examples.
One of the limitations is the fact that the computational cost of our method is proportional to the numberof basis functions. The level of detail in editing the flow
depends on the number of basis functions used. The
user can design a detailed flow field at the expense ofincreased computational and storage costs.
References
1. R. Bridson, Fluid Simulation for Computer Graphics,
Publisher, AK Peters, 2008.
2. R. Fattal, D. Lischinski, Target-driven smoke animation,
ACM Transactions on Graphics 23, 3, pp. 439-446, 2004.
3. Y. Kim, R. Machiraju, D. Thompson, Path-based controlof smoke simulaitons, In Proceedings of the 2006 ACM
SIGGRAPH/Eurographics symposium on Computer Animation, pp. 33-42, 2006.
4. M. B. Nielsen, B. B. Christensen, Improved variational
guiding of smoke aniemtions, Computer Graphics Forum29, 2, pp. 705-712, 2010.
5. A. Treuille, A. McNamara, Z. Popovic, J. Stam,
Keyframe control of smoke simulations, ACM Transactions on Graphics 22, 3, pp. 716-723, 2003.
6. T. D. Witt, C. Lessig, E. Fiume, Fluid simulation using
Laplacian eigenfunctions, ACM Transactions on Graphics
31, 1, Article 10, 2012.
7. Z. Yuan, F. Chen, Y. Zhao, Pattern-guided smoke animation with lagrangian coherent structure, ACM Transactions on Graphics 30, 6, Article 136, 2011.