Primal-Dual Optimization for Fluids

by   Tiffany Inglis, et al.

We apply a novel optimization scheme from the image processing and machine learning areas, a fast Primal-Dual method, to achieve controllable and realistic fluid simulations. While our method is generally applicable to many problems in fluid simulations, we focus on the two topics of fluid guiding and separating solid-wall boundary conditions. Each problem is posed as an optimization problem and solved using our method, which contains acceleration schemes tailored to each problem. In fluid guiding, we are interested in partially guiding fluid motion to exert control while preserving fluid characteristics. With our method, we achieve explicit control over both large-scale motions and small-scale details which is valuable for many applications, such as level-of-detail adjustment (after running the coarse simulation), spatially varying guiding strength, domain modification, and resimulation with different fluid parameters. For the separating solid-wall boundary conditions problem, our method effectively eliminates unrealistic artifacts of fluid crawling up solid walls and sticking to ceilings, requiring few changes to existing implementations. We demonstrate the fast convergence of our Primal-Dual method with a variety of test cases for both model problems.


page 1

page 6

page 7

page 9

page 10

page 13


A non-iterative domain decomposition method for the interaction between a fluid and a thick structure

This work focuses on the development and analysis of a partitioned numer...

A monolithic and a partitioned Reduced Basis Method for Fluid-Structure Interaction problems

The aim of this work is to present a brief report concerning the various...

Insights into the performance of loosely-coupled FSI schemes based on Robin boundary conditions

Robin boundary conditions are a natural consequence of employing Nitsche...

Optimization with nonstationary, nonlinear monolithic fluid-structure interaction

Within this work, we consider optimization settings for nonlinear, nonst...

Implementation, demonstration and validation of a user-defined wall-function for direct precipitation fouling in ANSYS Fluent

In a previous paper (Johnsen et al., 2015) and presentation (Johnsen et ...

Tracking via Motion Estimation with Physically Motivated Inter-Region Constraints

In this paper, we propose a method for tracking structures (e.g., ventri...

1 Introduction

Advances in fluid simulation have had a tremendous effect in engineering and graphics. Since fluids play an important role in our everyday lives, smoke and liquid simulations now routinely appear as regular elements in feature films, television and commercials. While engineering applications are primarily concerned with accuracy, the focus in graphics lies in simulating realistic behavior that can be controlled to tell a story. Therefore, it is important to provide controllable and visually plausible fluid solvers. Visual plausibility is crucial since people easily recognize unrealistic behavior due to their daily interactions with a wide range of fluid phenomena, e.g., pouring a glass of water, boiling a kettle, or driving past a chimney.

However, it is a recurring challenge to simultaneously achieve controllable and realistic behavior. Fluids are usually chaotic at human scales; miniscule perturbations regularly trigger large-scale behavior. It is almost impossible for artists to add small-scale details to a low-resolution simulation without changing the large-scale motion. Even more challenging is to modify the domain or rerun a simulation with different fluid parameters without affecting the large-scale motion. These problems can be mitigated by guiding the velocity within an optimization framework. We realize guiding by constraining the velocity at each time step to be arbitrarily close to the current and the target velocity. At the same time, we ensure that the resulting motion is divergence-free, and we allow for spatially varying guiding parameters. This framework also benefits animators in the special effects industry who may wish to create fictitious but visually plausible fluid motion to increase entertainment value.

Common fluid solvers typically feature boundary conditions (BCs) that lead to fluid unnaturally crawling up walls and sticking to ceilings, diminishing its visual plausibility. The culprit is the solid-wall BC that enforces a normal velocity of zero at obstacle walls, thus preventing fluids from separating. While it is physically correct to leave a thin film of fluid on the wall, its thickness in simulations is on the order of the discretization and usually far too big for realistic animations. One proposed remedy by Batty et al. [BBB07]

is to integrate inequality constraints into the pressure solve to allow for positive normal velocities at solid walls. Unfortunately, this greatly increases the complexity of the pressure solve of a fluid solver, and typically requires Quadratic Programming solvers. Our goal is to implement flexible separating BCs, while reusing the popular preconditioned conjugate gradient method. First, we classify boundary cells into separating and non-separating cells. We then enforce solid-wall BCs for the velocity at non-separating cells while ensuring zero divergence.

To guide simulations and realize separating BCs, we take inspiration from convex optimization. Gregson et al. [GITH14]

already established a connection between fluid pressure correction and convex optimization. This connection allows us to impose all required physical constraints on the velocity for both guiding and separating BCs via an efficient alternating minimization algorithm that exploits problem structure, removing the need for a monolithic and slow framework. We introduce a novel optimization scheme from the computer vision area, the

fast, or first-order Primal-Dual (PD) method proposed by Pock et al. [PCBC09]. It features an optimal convergence rate for non-smooth convex problems. Note that a whole class of primal-dual methods exists, but we will hereafter use the abbreviation PD to refer to this particular instance, which we will focus on due to its fast convergence properties.

Specifically, our contributions are

  • the introduction of a modular convex optimization approach for fluids based on the PD method;

  • a general method for handling the difficult problem of fluid guiding that involves spatially varying operators;

  • a fast method to approximately invert the linear system for flow guiding in arbitrary domains; and

  • a novel, practical way to handle separating BCs for liquids.

2 Previous Work

Fluid simulation has a long history in computer graphics [KM90]. One of the most widely used methods is the stable fluid solver [Sta99], for which many extensions have been proposed over the years, e.g., the grid-particle hybrid FLIP (Fluid Implicit Particle) method [ZB05], which is particularly popular for liquid simulations. A good overview of these methods is given in Bridson’s book [Bri08]. A key component of incompressible fluid simulation is ensuring that the time-evolved velocity is divergence-free (i.e., mass-preserving), which is typically implemented via Chorin-like pressure-projections, e.g.  [FF01] or  [BBB07] for variational approaches.

Convex Optimization

has become a powerful component of many computer vision algorithms. Even before, the works of Boyd et al. [BPC11] have laid the foundations for many popular optimization algorithms. Closer to computer graphics, Heide et al. [HRH13] have proposed using the PD method to correct the deficiencies of simple lenses. Recently, Heide et al. [HDN16] have solved the difficult task of choosing the optimal image priors and optimization algorithms for image processing tasks such as deconvolution, denoising or inpainting by applying PD as well. Iterated Orthogonal Projection (IOP), the algorithm proposed by Molemaker et al. [MCPN08], is a specialized method in the convex optimization area requiring all operators to be orthogonal projections. Approaches such as position-based fluids (PBF) [MM13]

also can be seen as an iterative application of a set of constraint projections. As the constraints are directly applied to the degrees of freedom, PBF is more closely related to IOP than PD, which achieves improved convergence by projecting the convex conjugate.

We extend upon the ideas presented by Gregson et al. [GITH14], where the method Alternating Direction Methods of Multipliers (ADMM) is applied, by introducing a more efficient splitting algorithm for flow guiding and BCs. Very recently, Narain et al. [NOB16] animated deformable objects with an ADMM algorithm combining a fast and robust nonlinear elasticity model with hard constraints. Our method can be seen as preconditioned version of ADMM, and its convergence is accelerated by an improved update direction for each iteration. O’Connor et al. [OV14] found ADMM outperforming PD for few iteration counts while using a restricted version of PD which does not exploit PD’s full potential of optimal update directions. We demonstrate the advantages of our approach over both ADMM and IOP in Section 4.1.


A limitation of fluid simulation for computer graphics is the issue of control. Simply changing the resolution of a simulation can alter its behavior significantly. A popular class of methods circumvents this problem by synthesizing smaller scales in a decoupled fashion [BHN07, KTJG08, PTSG09]. As a consequence, the details are easy to fine-tune, but do not tightly integrate with the base flow, which is left unmodified.

Fluid guiding is a challenging example of an inverse problem for fluids. Adjoint methods have been used in engineering and graphics [GP00, MTPS04] but require differentiation of the entire fluid solver. Guide shapes [NB11] are able to add detail to free surface flows by applying velocity conditions to high-resolution simulations distant to the interface but are limited in volumetric contexts. Other approaches aiming for high-level control have proposed the use of Lagrangian coherent structures [YCZ11] or sketches [PHT13] to give users intuitive controls. As our guiding technique takes an arbitrary flow field as input to calculate plausible and tightly coupled detail, such approaches would be a good complement for our method.

While Gregson et al. [GITH14]

also demonstrate preliminary results for guided flows, their approach employed the fast Fourier transform for filtering low-frequencies and pressure-projection. This is a very efficient approach, but it becomes impractical for more complex geometries and spatially varying guiding.

A simpler approach to guiding is via control forces [SY05b, FL04] or velocities [SY05a]. These approaches can result in artificially viscous flows, as noted by Thuerey et al. [TKRP06], who proposed a multi-scale approach based on control particles that controls low-spatial frequencies. Possibly the most similar approach to ours is the one by Nielsen et al. [NCZ09, NC10] who use a multi-scale formulation based on low pass filters. They arrive at a monolithic system with four degrees of freedom per cell, requiring a specialized solver to reduce run time. In contrast, our control scheme decouples into separate and more easily solved subproblems via recent developments in non-smooth optimization. Furthermore, for practical purposes, we introduce an upres method into our workflow [KTJG08, YCZ11, HMK11, RLL13, HK13] to facilitate the separation of low-resolution and high-resolution guiding.

Boundary Conditions

play a crucial role in fluid simulation, and as such have received a significant amount of attention. Foster et al. [FF01] describe how to allow for tangential motions of liquids, and are the first to note problems with liquid unnaturally sticking to domain boundaries. Batty et al. [BBB07] propose to solve inequality constraints with Quadratic Programming, while Chentanez et al. [CMF12] observe that they can incorporate these inequalities into their multigrid solver. Methods such as the FFT-based one of R. Henderson [Hen12] are similar in spirit to our method, as they separate the BCs from the divergence-free projection step, however, without targeting separating boundaries. Other approaches realize unilateral incompressibility to allow for separation effects [NGL10],[NGCL09], [GB13] but they also require complex solvers to solve their proposed Quadratic Programming Problems. We will demonstrate how to solve for separating motions with a regular conjugate gradient (CG) solver based on our modular optimization framework.

3 Methodology

In graphics, fluids are typically simulated by solving the incompressible Euler equations, written as


where is the flow velocity, the pressure, and the external body forces. Simulations commonly proceed via operator splitting to satisfy both constraints. First, using all but the pressure term in Eq. (1), an intermediate velocity field is computed. Then a pressure projection is applied to satisfy the divergence-free condition in Eq. (2). A pressure field that exactly counteracts the divergence is computed to correct the velocity field. Orthogonality of the curl-free and divergence-free components ensures that divergence-free components of the flow are not affected and allows the pressure projection to be interpreted as an Euclidean projection onto the space of divergence-free velocity fields.

The splitting approach closely resembles convex optimization approaches originally developed for imaging inverse problems and machine learning involving non-smooth or constrained objective functions. We show how several of these approaches can be adapted to solve difficult problems in fluids.

Convex Optimization

aims to solve problems of the form


where is a convex function that may be non-smooth or even discontinuous. If is the sum of two simpler functions (say ), then we have


A number of recently developed algorithms target this type of problem by employing an iterative divide-and-conquer approach. These algorithms are known as proximal methods and are defined in terms of so-called proximal operators (we use exclusively to denote the generic argument variable):


One such algorithm is the PD method [PCBC09], which solves a slightly more general problem:


for some linear operator . PD solves the problem iteratively by providing a series of variable updates that terminate when converges to the solution. The combination of and ensures that converges to the optimal value of Eq. 6. The updates are given by


where are parameters that affect convergence, and are the convex conjugates of and , respectively. For our problem, is simply the identity, which leads to . Note that if additionally , the iterative update scheme reduces to ADMM. A more appropriate choice () leads to optimal control over convergence. As for , it is not necessary to compute it directly. The proximal operator can be transformed using Moreau’s identity:


The variable updates are thus reduced to


A more in-depth discussion of PD can be found in [CP11].

The advantage of using proximal methods is that the optimization can be performed separately for the two objective functions, allowing difficult optimizations to be split into more manageable components. Also, depending on the form of and , many special cases can be significantly simplified [BPC11] by exploiting their mathematical structures.

To the best of our knowledge, PD has not yet been applied to fluid problems, despite its provably optimal convergence properties for the class of problems of Eq. (6). A pseudocode implementation of our PD-based optimization method is given in Algorithm (1

) for one time step. Comparing to previous methods, the improved convergence is achieved with only a minimal increase in computational cost. A few more vector additions and multiplications are necessary, which typically are negligible compared to the cost of the proximal operators. We demonstrate the convergence of our method and compare it to other methods in Sections 

4.1 and 5.1. For reference, we briefly review IOP and ADMM in Appendix B.

Convex Optimization of Fluids

In fluid simulation, we often encounter problems of the form

subject to

where is the velocity field we seek, and is the space of divergence-free velocity fields. must be a convex function. Practically, this can be either a quantity we are trying to minimize, or a hard constraint that must be satisfied (in which case, would be an indicator function).

The second constraint requires to be divergence-free. Removing the divergent part of the flow can also be viewed as an orthogonal projection [Cho68, PB13]. Gregson et al. [GITH14] made the key observation that the proximal operator for can be easily computed via a pressure projection. In other words, we have


where denotes a projection onto with a commonly used Poisson solver. Hence, formulating a fluid problem this way allows the optimization algorithm to be easily integrated into a common fluid solver—we simply replace the call for a pressure projection subroutine with a call to the PD optimization step outlined in Algorithm (1). We check for convergence using a threshold parameter , and stop the algorithm once the per-iteration change of falls below this threshold (Algorithm (1), line 12). Note that we define the pressure projection to include only the calculation of the pressure values, and a subtraction of the pressure gradient from , excluding any optional modifications of the velocity field.

1:procedure PD(, , , , )
2:     while  do
3:         // -update
5:         // -update (using adaptive CG accuracy)
7:         // -update
9:         // check stopping criterion
12:         if then break      
13:     return
Algorithm 1 Our PD-based method for fluid simulation

The pressure projection implemented by a CG solver is usually the most expensive part in regular fluid simulation, and its effect is magnified in our algorithm due to its iterative nature. However, this iterative set-up gives us an opportunity to apply an adaptive scheme. Let be the accuracy of the CG solver. We choose its value starting with a high threshold (e.g., ) and then adaptively decrease it over the PD iterations. We decrease as soon as the per-iteration change of is close to ; until reaches the desired final accuracy (e.g., ). Using this adaptive CG scheme greatly accelerates the performance by cutting down on CG iterations in the beginning of the optimization when the divergence-free constraint does not need to be strictly enforced.

Algorithm (1) summarizes the general framework of our method as it applies to fluid simulation. Specific applications call for different definitions of , which in turn affects how the -update is computed. In the following sections, we discuss how to apply this method to the fluid guiding and the separating BC problem. For each application, we define the appropriate and discuss how to compute its proximal operator.

4 Fluid Guiding

In fluid guiding, the goal is to minimize the change applied to the current velocity field such that the resulting velocity field follows the large-scale motions of a given target velocity. The objective function is given by


where is the current velocity field (after advection and before pressure projection), is the target velocity field, is the guided velocity field, is the guiding weights matrix, and is the Gaussian blur matrix. The first term of this objective function is minimal for a solution that matches the target velocity when blurred by , while the second term penalizes solutions far away from the current flow field.

In order to keep the application general, both matrices in our objective function are spatially varying (Nielsen and Christensen [NC10] also performed fluid guiding with spatially varying guiding weights, but not with spatially varying blur). is a diagonal matrix containing spatially varying weights that control the guiding strength (larger entries denote weaker guiding), and is the Gaussian blur matrix that applies a blur of radius only to fluid cells so that we can handle boundaries and obstacles in our domain.

Our objective function in Eq. (16) is quadratic. That is, it can be expressed as




Note that since is symmetric.

We can combine (three equations per cell in 3D) and the divergence-free constraint (one equation per cell) into a linear system . But since this system is overconstrained, we solve it in the least-squares sense by considering


where and . This quadratic system can be solved using an iterative solver such as the CG solver. However, this approach is infeasible since the number of elements in the matrix grows by , where is the volume of the simulation domain. Thus PD is still a good choice even for seemingly simple quadratic energies. We will see later how this direct CG solver compares to our algorithm in terms of performance.

Instead, we make use of the quadratic property of to simplify its proximal operator to


Factoring out , we obtain



(a) Target velocity
(b) Linear velocity blend
(c) Detail-preserving blend
Figure 1: Naïve guided simulations with circular targets.

Computing for every iteration is slow. Instead, we precompute and since they do not rely on previous iterations. In fact, using the Sherman-Morrison-Woodbury formula, can be approximated as


See Appendix C for derivation details.

Lastly, spatially invariant Gaussian blurs are symmetric, so . They are also separable, and can be applied independently in each dimension. This combined with the matrix inversion approximation greatly speeds up the -update. For Gaussian blurs with spatial variation, we still use a symmetric and hence separable Gaussian blur (with different blur radii) for each cell, but itself is no longer symmetric. However, if the spatial variation is limited (e.g. different blur radii on the left and right sides of the simulation domain), we approximate by , since for regions of constant blur. We found our approximation to work well in practice, and we will show examples using this approximation later. Our PD guiding scheme is independent of the particular choice and implementation of the blur kernel. Given a fast implementation for calculating , it would be straight-forward to extend our method with a full evaluation of .

Additionally, our approach for applying the blur kernel easily allows for interior boundaries by using a blur radius of zero for obstacle cells. This is a key difference from guiding scheme based on the fast Fourier transform [GITH14], which typically require periodic domains without internal boundaries.

Algorithm (2) summarizes how is approximated. Note that the two extra parameters and in Algorithm (2) as compared to Algorithm (1) are specific to the guiding application.

1:procedure proxF(, , , )
2:      // can be precomputed
3:      // inverse of diagonal matrix
4:     return
Algorithm 2 Approximation for in fluid guiding


In this section, we will first demonstrate the efficacy of our method using 2D examples, followed by more impressive 3D results. To simplify the notation, when the guiding weight is spatially invariant (say, everywhere), we will write . And when the guiding weight is spatially varying (say, and on the left and right side on the simulation domain, respectively), we will write . Similarly, spatially varying blur radii would be written as .

We first compare our method to two naïve guiding methods, with a counterclockwise circular velocity field as the target , as shown in Figure 0(a). Linear velocity blend computes the new velocity as a linear combination of the current velocity and the target velocity: , where . Detail-preserving guiding [TKRP06] subtracts the large-scale (i.e., blurred) motions from the current velocity before adding the target velocity: . The idea is to create a new velocity with large-scale motions from and small-scale details from . Both methods are followed by a pressure solve to ensure zero divergence.

With linear velocity blend, guiding strength is increased with smaller . The method’s main weakness is that small fluid details are smoothed out with strong guiding, and more details come at the expense of reduced trajectory control (see Figure 0(b)).

Detail-preserving blend allows guiding of large-scale motions towards the target velocity while preserving details. However, the details are difficult to control and have an unnatural frosted glass appearance (see Figure 0(c)).

In contrast, our method has much more flexible motion control when applied to fluid guiding, with affecting the large-scale guiding strength and the blur radius controlling the small-scale details. Figure 2 shows how the parameters affect guiding for a 2D smoke simulation with a circular target velocity. Notice that larger values allow for more freedom to deviate from the target, while larger values lead to the formation of larger vortices.

We also compare our algorithm to wavelet turbulence [KTJG08], as shown in Figure 3. Although wavelet turbulence is fast—since it is purely a postprocessing technique—the resulting vortices do not couple as tightly and realistically as with our method.

Figure 2: Guided simulations with circular targets, using our method with varying and for multi-scale motion control.
Figure 3: Upsampling from to using wavelet turbulence (top) and using our method (bottom). Our method yields significantly more coherent vortices.


Next, we examine our method in relation to other proximal methods, namely, ADMM and IOP (see Appendix B for details). ADMM is fairly straightforward to apply due to its similarity to PD. IOP can be interpreted as a fluid-specific version of the Projection onto Convex Sets (PoCS) algorithm [BB96], which solves convex feasibility problems (locating an intersection point of convex sets) rather than the more general problem solved by ADMM and PD. Applying IOP naïvely alternates between the minimizer of and its closest divergence-free neighbor but does not yield the divergence-free minimizer of , as illustrated by the failure case in Figure 4. In practice, we discovered instances for which IOP found plausible approximate solutions in a short time, but its reliability is limited for general guiding applications. As such, we focus on only comparing our method to ADMM. All performance measurements were done on PCs with Intel Xeon E5-1650 CPUs.

Figure 4: A divergent target velocity leads to IOP failure, while ADMM and PD converge successfully. In this case, even though the IOP solution is divergence-free, it lands on an isocontour farther from the center than the PD/ADMM solution, which is the true minimum.

For fairness of comparison, we experimented with several test scenes to deduce as-optimal-as-possible parameter sets for both ADMM and our method. These parameters, on average, produce the fastest convergence rates under various guiding weights and blur radii. The optimal parameters are correlated with the guiding weight; we analyzed their relationship and define the parameters in terms of the average guiding weight, , which is simply the mean of the diagonal terms of . For our method, we chose , and for ADMM, . Before comparing their performance, it is important to note that we apply one of our contributions crucial for performance (the matrix inverse approximation and the separable Gaussians) to both methods.

(a) 2D circular example
(b) 3D obstacle example
Figure 5: Mean number of iterations per time step for guiding.

We compare the performance of ADMM and our method for 2D and 3D problems. The 2D problem is a simulation guided by a circular target velocity, similar to the one shown in Figure 2, but with spatially varying weights and . Specifically, the guiding weight is fixed at 1 for the right side of the domain, while the left side takes values from . Figure 4(a) compares, for the two methods, the mean number of iterations required to reach convergence. When the spatial variation is low, both methods perform decently. However, as the spatial variation increases, ADMM takes much longer to converge.

The 3D problem is a simulation with an obstacle, similar to the one in Figure 12, but with the same and parameters as the 2D problem. The performance results are shown in Figure 4(b). Again, our method outperforms ADMM. In fact, the difference in convergence rates is even more apparent for the 3D example for large spatial variations. The run times for both ADMM and our method very closely correlate with the mean iteration counts. For brevity, the run times are given in Table 1.

Guiding weight W
(2,1) (4,1) (8,1) (16,1)
2D Our method 0.14 0.17 0.26 0.38
ADMM 0.26 0.41 0.76 1.89
3D Our method 16.6 19.5 26.3 38.4
ADMM 30.4 67.8 176.9 325.6
Table 1: Mean run time (in seconds) per time step for guiding.
Figure 6: Performance of our method on the 3D tornado example at various resolutions, including (a) the mean run time per time step, and (b) the mean run time per time step per grid cell.

In addition to the comparison with ADMM, we also analyze how our method scales as the resolution increases. In Figure 5(a), we run the guided 3D tornado simulation shown in Figure Primal-Dual Optimization for Fluids at five different resolutions and compare the mean run time per time step. The scaling factor is with respect to the lowest resolution ; for instance, a scaling factor of 3 corresponds to the resolution . Notice that even at the largest resolution , each time step takes less than four minutes to complete. However, since the run times appear to scale exponentially with respect to resolution, we plot the mean run time per grid cell (see Figure 5(b)) to show that the scaling is in fact linear.

Previously, we mentioned that applying a CG solver also works for the guiding problem, although its poor performance renders it infeasible. To test this, we ran a small example with a circular target velocity field and found the direct CG solver (see Eq. (4)) approach to be 4000 times slower than our method.

The performance tests above are all done with our generic guiding method, which handles spatially varying weights and blurs. In the special case with spatially invariant weights, we can apply a non-iterative method to greatly speed up the simulation. The method involves defining in terms of the non-divergent components of the input velocities and then finding its minimizer. Since differentiation commutes with convolution, the result will be automatically divergence-free for spatially invariant weights. Although this non-iterative method offers potential speed-up, we focus on the general case of spatially varying operators that requires iterating.

4.1 Guiding Results

(a) with up-res
(b) without up-res
Figure 7: Simulations guided by a star-shaped target.

To realize interesting guided flows, we found an iterative up-res workflow that works well in practice. This is in line with previous work, where a low-resolution input is refined in subsequent stages to yield a final result [KTJG08, YCZ11, HMK11, RLL13, HK13]. Examples of this process are shown in Figure 6(a). We start with a simulation guided by a synthetic star-shaped velocity field with and . Then we run high-resolution simulations with various values. If the same guiding is done without up-res, that is, if we run the guided simulation directly as shown in Figure 6(b), then we lose the ability to guide the fluid to desirable shapes at multiple levels.

We further demonstrate this up-res process for 3D simulations. Figure 10 shows a simple plume example. We first run a low-resolution () simulation to capture the velocities. Then a resimulation is performed, guided by the upsampled velocity field. The blur radii can be kept spatially invariant or spatially varying for different effects. In particularly, increasing the blur radius introduces more small-scale details. In our spatially varying example, we use a blur kernel with a sharp transition between two blur radii to contrast between the two guided regions. If desired, this obvious seam can be softened via a more gradual transition between the two blur radii. Note that even though the simulation uses the approximation , it works quite well in practice.

Our second example in Figure 10 shows 3D simulations following star-shaped input velocities. We use a target velocity as guide for a final resolution, with . Here we make use of spatially varying weights, with on the left side of the domain, and on the right. The precomputation time is negligible compared to the overall run time: around seconds per time step for a simulation.

Next, we show a tornado simulation. We first run a low-resolution () simulation using a cylindrical target velocity with a small upward component. Once the general shape is fixed, a resimulation is performed, guided by the upsampled velocity field. The results in Figure 10 show the effectiveness of varying to achieve different levels of turbulence.

Figure 9: A star-shaped guided simulation upsampled from a guided simulation with spatially varying weights.
Figure 8: A simple plume simulation (left) is upsampled to with constant blur (middle, ) and spatially varying blur (right, ).
Figure 9: A star-shaped guided simulation upsampled from a guided simulation with spatially varying weights.
Figure 10: A guided simulation upsampled from a guided simulation with varying .
Figure 8: A simple plume simulation (left) is upsampled to with constant blur (middle, ) and spatially varying blur (right, ).
Figure 11: A simulation with an obstacle (inset) is upsampled to with and .

Finally, Figure 11 demonstrates an example with obstacle, for which we upsampled from to with and . This high-resolution version captures the input motion, but develops many interesting small-scale details. Figure 12 shows a similar simulation with and spatially varying weights, where the left side has while the right side has . The result has little detail on the left side due to strong guiding from the low-resolution simulation, whereas the right side behaves like a regular smoke simulation due to weak guiding. This guided simulation with an obstacle is a case that cannot be handled by previous Fourier-domain guiding schemes [GITH14] due to the interior boundaries. Arbitrary boundaries are easily handled by our separable approximation to the fluid guiding proximal operator while remaining highly efficient and easily parallelizable.

Figure 12: Guided obstacle scene with spatially varying weights.

5 Separating Solid-Wall BCs

Another application of our PD-based method is the realization of separating solid-wall BCs. As a motivational example, consider Figure 13, which shows a 2D breaking dam simulation. When generated by a common CG pressure solver (top row), the fluid exhibits the undesirable behavior of sticking to the ceiling and getting stuck in corners. In contrast, our method (bottom row) features a clean separation of the fluid from all solids. In many visual effect settings that aim for large-scale fluid flow, the separation is preferable due to its improved realism.

We achieve this behavior by implementing separating solid-wall BCs in the form of an inequality constraint (with and denoting velocity and obstacle normal) without transforming the linear system of equations for the pressure solve into a more complicated problem. As outlined in Section 3, we split the problem into two simpler objective functions, with controlling the BCs and enforcing zero divergence (see Eq. (15)). With this setup, a regular CG solver is employed to compute , while is handled by an efficient projection and classification scheme, as described next.

Velocity BCs

The proximal operator ensures that velocities at obstacle surfaces never point into the solid (i.e., negative normal velocity components at fluid-solid faces). Separating velocities (i.e., positive normal components), on the other hand, are allowed and thus left alone. We assume an obstacle is static and has a velocity of zero at its surface.

Figure 13: Breaking Dam, Non-Separating (top) versus Separating BCs (bottom) at .

BCs that form a linear constraint on the velocity field, like requiring velocity components to be zero, can be expressed as orthogonal projections [MCPN08]. Therefore, reduces to a projection onto the space of velocities without flows into obstacles. In its simplest form, this proximal operator sets normal velocity components to zero wherever . The pressure solver, on the other hand, is unaware of any obstacle boundaries and uses free surface Dirichlet BCs at liquid surfaces. Later, we present an accelerated version that makes the pressure solver partially aware of obstacles. For now, all obstacles are fully handled by .

Although this simple velocity projection works for exact solutions, it breaks down when solving up to finite accuracy. The pressure solver can accumulate small positive velocities during iterations of our optimization procedure, leading to small separating motions where the liquid should only be standing still or moving tangentially. This can happen for hydrostatic cases, for instance, where the final velocity at the wall should be exactly zero. We address this problem by introducing a classification step with hysteresis.

We assume by default that our cells are separating boundary cells with a Dirichlet condition, and then create a list of cells that are explicitly not allowed to separate. For the classification, in each iteration of our algorithm, we accumulate all motions into a wall returned by the CG solver for a cell and store it in . This allows for a solid cell classification with temporal coherence. A non-separating cell only becomes separating if the magnitude of its velocity component away from the wall exceeds the magnitude of . Additionally, a change in the separation classification is only made if the absolute value of is greater than the accuracy of the CG solver . Otherwise, the cell keeps its previous state. The last two steps effectively implement a hysteresis that prevents cells from changing status due to numerical errors, which is crucial for a reliable convergence of our method.

Figure 14 illustrates this process. It shows two possible fluid behaviors at the face of a solid: non-separating (left) and separating liquid (right). Both velocities fulfill (black bars), and can thus potentially change state. On the left of Figure 14, the cell face has a large accumulated value (red arrow). Assuming that the face velocity (green arrow) exceeds the threshold , the cell is classified as non-separating, and is later on set to zero despite its positive value. In contrast, a positive normal velocity component greater than indicates that fluid is currently separating from the solid cell. The liquid boundary is treated as a free surface in this case, and its velocity is left unmodified.

Figure 14: Classification of solid cells into non-separating (left) and separating (right). The green and red arrows denote and respectively, while the black bars indicate .

Algorithm Summary

Algorithm (3) shows pseudocode for the classification and the proximal operator for separating BCs. The classification proceeds as outlined above, updating the list of non-separating cells based on the velocity field . The normal velocity components at non-separating cell faces are set to zero in the proximal operator . The variable is a generic argument variable being a combination of multiple velocity fields, including dual variables.

The PD algorithm solving for separating BCs follows the generic PD implementation outlined in Algorithm (1). We exchange the generic proximal operator by our BC specific projection from Algorithm (3). Additionally, we call the classification function with after line 6 of Algorithm (1).

1:procedure classify(, , , )
2:     for all boundary faces , with  do
3:         if  then
4:                  // non-separating
6:         else if ( then
7:                  // make separating
9:procedure proxF(, )
10:     // velocity projection for marked cells
11:     for all do
12:     return
Algorithm 3 Classification and for BCs

We adopt two smaller modifications from previous work that improve convergence. The first is a Krylov method from IOP [MCPN08], which we identified to work nicely within our BC solver (the full algorithm can be found in Appendix B). Additionally, we use the adaptive parameter scheme [CP11] for dynamic values of , and . The parameter updates are the following after choosing a suitable value for , and (we use , , ): , , and .

Both methods reduce the overall run time, but are specific extensions for proximal operators given by orthogonal projections. As such, we use the extensions for all BC problems, but not for our guided simulations.

Accelerated BC Solver

So far, the solid-wall BC handling is fully done in while simply ensures incompressibility. This is in line with the standard splitting approach in other methods [MCPN08, Hen12], and we demonstrate in the next section that our optimization scheme very efficiently calculates the solution in this case.

However, in practice, the total number of iterations can be reduced significantly by letting the pressure solver take care of all BCs that it is capable of enforcing correctly: retaining the input normal velocity components at fluid-solid faces by enforcing a zero pressure gradient at walls. For this version of our solver, we update the BCs of the pressure solver in in accordance with our classification: Neumann BCs for non-separating cells, and Dirichlet (free surface) BCs for separating ones. This effectively locks the classification after few iterations, leading to a stable solution within the next iteration. However, it prevents non-separating cells from changing their state back to separating ones. Thus, the accelerated BC solver does not yield the same result as the standard version, but achieves much higher performance while featuring only negligible differences. A summary of our accelerated solver with pseudo-code can be found in Appendix D.

5.1 Evaluation and Results

We now evaluate the performance of our BC solver, and demonstrate the importance of allowing liquid to separate from walls in high-resolution simulations.111 Thresholds for our algorithm are set to (as liquids are more sensitive to mass loss) and . To validate that our PD solver yields the correct results, we use our method to calculate regular non-separating boundaries for a breaking dam setup, which can also be handled by a regular CG solver. In this case, we can achieve arbitrary accuracies depending on the choice of parameters for CG and our scheme. These experiments show that our method converges to the correct solution.

For separating boundaries, we can no longer compare to a regular CG solver. Instead, we evaluate the performance of our method against methods from previous convex optimization work: IOP and ADMM. IOP has been applied to enforce non-divergence and complex BCs simultaneously in [CMF12]. For the standard BC solver simulating a 2D breaking dam, our method converges six times faster than ADMM, and more than twice as fast as IOP, as shown in Figure 14(c). Figure 14(a) and Figure 14(b) show the mean number of IOP/ADMM/PD and CG iterations. Our method generally requires a lower number of PD iterations compared to ADMM and IOP. Due to our adaptive CG accuracy scheme, the number of PD iterations is increased while the total amount of CG iterations is decreased. Both numbers influence the performance, but the total number of CG iterations more strongly influences runtime as shown in Figure 14(c). Similar behavior is observed for a 3D complex breaking dam simulation (see Figure 17 for a visual example). The run time measurements are shown in Figure 16 where ADMM is omitted due to its impractical run time. The accelerated solver speeds up all methods significantly. In this case, our PD-based method performs on par with IOP, which we attribute to the low number of iterations required in this case; the higher-order convergence of the PD method does not develop its full potential in such a setting.

(a) IOP/ADMM/PD Iter
(b) Sum of CG Iter
(c) Runtime in seconds
Figure 15: 2D breaking dam, mean quantities per time step.
Figure 16: 3D complex breaking dam, mean run time per time step.

Figure 17 highlights that regular non-separating walls often yield undesirable results in practical settings. We simulate a complex 3D scene with liquid splashing onto multiple obstacles with a resolution of . A regular solver leads to large amounts of liquid crawling along ceilings, and sticking to walls, while our separating boundaries give a much more believable large scale look, with liquid naturally separating from obstacle boundaries. For this setting, our accelerated BC solver requires on average per frame. Comparing the overall time to generate this result (including surface generation) to a version with a regular pressure solver, our method increases run time by only , with the added benefit of enabling separating BCs. This is a very practical result, considering that it was achieved based on a regular CG solver, without the need for specialized methods, such as Linear Complementarity Problems solvers [GB13]. Although solving LCPs has been studied in detail, these methods still fall into complexity classes that make large scale solves infeasible. An indication for this can be found in [GB13], where solving times of ca. 25s per LCP solve were given for a example.

Figure 17: A scene with complex obstacles. Regular, non- BCs lead to the liquid sticking to walls (left), while our solver allows for naturally separating liquids (right).

6 Discussion and Conclusion

We presented a general framework for incorporating the Primal-Dual method into fluid simulation, and demonstrated two applications: fluid guiding and separating BCs. We proposed a generic version of the Primal-Dual based optimization scheme with fast convergence for general proximal operator subproblems. In addition, we discussed several extensions that are particularly well-suited for optimizations with orthogonal projections as subproblems. Additionally, we demonstrated a novel formulation for the flow guiding problem, and an efficient approach for simulating liquids with separating BCs.


One limitation of our fluid guiding method is a lack of shape controls. Unlike smoke, liquids can require shape constraints to achieve a desired outline or shape. We have focused on velocity guiding in this work, so controlling the shape of liquids will require extensions that take the position of the liquid’s surface into account. Such constraints should integrate well into our overall pipeline, and we plan to investigate this topic as future work. Another area of improvement is bridging the gap between our optimization framework and an artist’s workflow. Although we gave examples of different target velocities used for guiding, there is still a lot more to explore in terms of how to achieve various artistic visual effects intuitively. A limitation of our accelerated separating BC solver is that it can lead to slight deviations from the accurate, standard solution. We have not encountered visual artifacts resulting from these inaccuracies, and we believe that the improved performance typically outweighs these slight deviations. The specialized multigrid solvers [CMF12] could possibly outperform our BC solver. We believe that the attractiveness of our BC solver stems from its modularity and fast convergence. Adding support for wall separating BCs given an existing pressure solver requires little code with our method.


Our method is a very generic approach applicable to a large range of problems in fluid simulation. The price of this generality is that it may not be as fast as specialized methods tailored to specific problems. However, the modularity of our approach makes it easy to incorporate into existing implementations. As such, it has the potential to add powerful functionality, such as high-level flow guiding, into existing solvers without the need for complicated extensions.

Furthermore, there is a large number of interesting avenues to be explored with high-level optimizations of fluids flows. For example, we are interested in exploring shape optimization to adapt the geometry of obstacles with respect to their flow properties, and partial resimulations of regions in a flow.

7 Acknowledgments

This work was funded by the ERC Starting Grant realFlow (StG-2015-637014) and the NSERC Postdoctoral Fellowships Program.

Appendix A Notation

For reference we provide a quick description of our notation. More detailed descriptions can be found in the body of the paper.

  • : first objective function in PD (application-dependent)

  • : second objective function in PD for incompressibility

  • : general input velocity field for several algorithms, a combination of several velocity fields, possibly dual variables

  • , , : variables with iterative updates in PD

  • , , : PD parameters controlling convergence rate

  • , , , , , , , : variables specific to fluid guiding

  • : fluid velocity for the BC problem

  • : normal of a solid cell

  • : memory velocity field

  • : set of non-separating fluid-solid faces

Appendix B ADMM and IOP

Iterated Orthogonal Projection (IOP) [MCPN08]—a method similar to von Neumann’s alternating projections [BPC11]—requires both subproblems to be expressed as orthogonal projections


The Krylov method can improve its convergence rate:

1:procedure krylov(, , , )
3:     if  then
4:           // correction vector
8:          if then      

The Alternating Direction Method of Multipliers (ADMM) [BPC11, GOSB14] is a proximal method more general than IOP given by


Instead of the three parameters in PD that control convergence rate, ADMM only has one such parameter, .

Appendix C Inverse Matrix Approximation for Fluid Guiding

Here we present the details for deriving the inverse matrix approximation in Eq. (26).

We want to invert a matrix of the form , where contains the large diagonal terms and the small off-diagonal terms. By the Sherman-Morrison-Woodbury Formula,


Now since and both contains small value entries, is approximately zero. Hence


Note that the calculation of is trivial since is diagonal.

Appendix D Accelerated BC Solver

Below, we summarize our accelerated solver for separating solid-wall BCs. Instead of fully handling the solid-wall BCs in a proximal operator outside of the pressure projection, our accelerated solver employs the commonly used Neumann BCs during the pressure projection depending on the current classification of solid-wall cells. First, all solid wall cells are classified as separating (i.e. Dirichlet BCs for the pressure solver). Depending on the initial velocity field, solid wall cells with are classified as non-separating cells (Neumann BCs for the pressure solver). We then iteratively set the velocity BCs, apply a regular CG pressure solver and re-classify solid cells with the classify procedure.

In our standard BC solver (Section 5), a memory field is used to determine whether a cell is allowed to change its state to separating. In the accelerated case, is not used, instead the pressure solver enforces solid wall BCs for non-separating cells. This prevents the normal velocity from becoming positive, and from changing its state back to separating. In the following procedure, is replaced with the placeholder "".

1:procedure solvePressureWithSeparatingBCs(, )
2:      // no cell is classified as non-separating yet
3:     classify(, , , )
4:     while  did not change in last call of classify do
5:          (, )
6:          (, )
7:          classify(, , , )      


  • [BB96] Bauschke H. H., Borwein J. M.: On projection algorithms for solving convex feasibility problems. SIAM review 38, 3 (1996), 367–426.
  • [BBB07] Batty C., Bertails F., Bridson R.: A fast variational framework for accurate solid-fluid coupling. ACM Trans. Graph. 26, 3 (July 2007).
  • [BHN07] Bridson R., Houriham J., Nordenstam M.: Curl-noise for procedural fluid flow. ACM Transactions on Graphics (TOG) 26, 3 (2007), 46.
  • [BPC11] Boyd S., Parikh N., Chu E., Peleato B., Eckstein J.: Distributed optimization and statistical learning via the alternating direction method of multipliers. Found. Trends Mach. Learn. 3, 1 (Jan. 2011), 1–122.
  • [Bri08] Bridson R.: Fluid Simulation for Computer Graphics. AK Peters/CRC Press, 2008.
  • [Cho68] Chorin A. J.: Numerical solution of the Navier-Stokes equations. Mathematics of computation 22, 104 (1968), 745–762.
  • [CMF12] Chentanez N., Mueller-Fischer M.: A multigrid fluid pressure solver handling separating solid boundary conditions. Visualization and Computer Graphics, IEEE Transactions on 18, 8 (2012), 1191–1201.
  • [CP11] Chambolle A., Pock T.: A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision 40, 1 (2011), 120–145.
  • [FF01] Foster N., Fedkiw R.: Practical animation of liquids. In Proceedings of ACM SIGGRAPH (2001), ACM, pp. 23–30.
  • [FL04] Fattal R., Lischinski D.: Target-driven smoke animation. In ACM Transactions on Graphics (TOG) (2004), vol. 23, ACM, pp. 441–448.
  • [GB13] Gerszewski D., Bargteil A. W.: Physics-based animation of large-scale splashing liquids. ACM Trans. Graph. 32, 6 (2013), 185.
  • [GITH14] Gregson J., Ihrke I., Thuerey N., Heidrich W.: From capture to simulation: connecting forward and inverse problems in fluids. ACM Trans. Graph. 33, 4 (2014), 1–11.
  • [GOSB14] Goldstein T., O’Donoghue B., Setzer S., Baraniuk R.: Fast alternating direction optimization methods. SIAM Journal on Imaging Sciences 7, 3 (2014), 1588–1623.
  • [GP00] Giles M. B., Pierce N. A.: An introduction to the adjoint approach to design. Flow, turbulence and combustion 65, 3-4 (2000), 393–415.
  • [HDN16] Heide F., Diamond S., Nießner M., Ragan-Kelley J., Heidrich W., Wetzstein G.: ProxImaL: Efficient Image Optimization using Proximal Algorithms. ACM Transactions on Graphics (TOG) 35, 4 (2016).
  • [Hen12] Henderson R. D.: Scalable fluid simulation in linear time on shared memory multiprocessors. In Proceedings of the Digital Production Symposium (2012), DigiPro ’12, ACM, pp. 43–52.
  • [HK13] Huang R., Keyser J.: Automated sampling and control of gaseous simulations. The Visual Computer 29, 6 (2013), 751–760.
  • [HMK11] Huang R., Melek Z., Keyser J.: Preview-based sampling for controlling gaseous simulations. In Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (2011), SCA ’11, ACM, pp. 177–186.
  • [HRH13] Heide F., Rouf M., Hullin M. B., Labitzke B., Heidrich W., Kolb A.: High-quality computational imaging through simple lenses. ACM Trans. Graph. 32, 5 (Oct. 2013), 149:1–149:14.
  • [KM90] Kass M., Miller G.: Rapid, stable fluid dynamics for computer graphics. In SIGGRAPH (1990), ACM Press, pp. 49–57.
  • [KTJG08] Kim T., Thuerey N., James D., Gross M.: Wavelet Turbulence for Fluid Simulation. ACM Trans. Graph. 27, 3 (2008), article 50.
  • [MCPN08] Molemaker J., Cohen J. M., Patel S., Noh J.: Low viscosity flow simulations for animation. In Proceedings of the Symposium on Computer Animation (2008), Eurographics Association, pp. 9–18.
  • [MM13] Macklin M., Müller M.: Position based fluids. ACM Transactions on Graphics (TOG) 32, 4 (2013), 104.
  • [MTPS04] McNamara A., Treuille A., Popović Z., Stam J.: Fluid Control Using the Adjoint Method. ACM Trans. Graph. 23, 3 (2004), 449–456.
  • [NB11] Nielsen M. B., Bridson R.: Guide shapes for high resolution naturalistic liquid simulation. ACM Transactions on Graphics (TOG) 30, 4 (2011), 83.
  • [NC10] Nielsen M. B., Christensen B. B.: Improved variational guiding of smoke animations. In Computer Graphics Forum (2010), vol. 29, Wiley Online Library, pp. 705–712.
  • [NCZ09] Nielsen M. B., Christensen B. B., Zafar N. B., Roble D., Museth K.: Guiding of Smoke Animations through Variational Coupling of Simulations at Different Resolutions. In Proceedings of the Symposium on Computer animation (2009), pp. 217–226.
  • [NGCL09] Narain R., Golas A., Curtis S., Lin M. C.: Aggregate dynamics for dense crowd simulation. In ACM Transactions on Graphics (TOG) (2009), vol. 28, ACM, p. 122.
  • [NGL10] Narain R., Golas A., Lin M. C.: Free-flowing granular materials with two-way solid coupling. ACM Transactions on Graphics (TOG) 29, 6 (2010), 173.
  • [NOB16] Narain R., Overby M., Brown G. E.: ADMM projective dynamics: Fast simulation of general constitutive models. In Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation (2016), Eurographics Association, pp. 21–28.
  • [OV14] O’Connor D., Vandenberghe L.: Primal-dual decomposition by operator splitting and applications to image deblurring. SIAM Journal on Imaging Sciences 7, 3 (2014), 1724–1754.
  • [PB13] Parikh N., Boyd S.: Proximal Algorithms. Foundations and Trends in Optimization 1, 3 (2013), 123–231.
  • [PCBC09] Pock T., Cremers D., Bischof H., Chambolle A.: An algorithm for minimizing the mumford-shah functional. In Computer Vision (2009), IEEE, pp. 1133–1140.
  • [PHT13] Pan Z., Huang J., Tong Y., Zheng C., Bao H.: Interactive localized liquid motion editing. ACM Trans. Graph. 32, 6 (2013), 184:1–184:10.
  • [PTSG09] Pfaff T., Thuerey N., Selle A., Gross M.: Synthetic turbulence using artificial boundary layers. In ACM Transactions on Graphics (TOG) (2009), vol. 28, ACM, p. 121.
  • [RLL13] Ren B., Li C.-F., Lin M. C., Kim T., Hu S.-M.: Flow field modulation. IEEE Transactions on Visualization and Computer Graphics 19, 10 (2013), 1708–1719.
  • [Sta99] Stam J.: Stable Fluids. In Proceedings ACM SIGGRAPH (1999), pp. 121–128.
  • [SY05a] Shi L., Yu Y.: Controllable smoke animation with guiding objects. ACM Trans. Graph. (TOG) 24, 1 (2005), 140–164.
  • [SY05b] Shi L., Yu Y.: Taming Liquids for Rapidly Changing Targets. In Proceedings of the ACM SIGGRAPH/Eurographics symposium on Computer animation (2005), pp. 229–236.
  • [TKRP06] Thuerey N., Keiser R., Ruede U., Pauly M.: Detail-Preserving Fluid Control. In Proceedings of the Symposium on Computer animation (2006), pp. 7–12.
  • [YCZ11] Yuan Z., Chen F., Zhao Y.: Pattern-guided smoke animation with lagrangian coherent structure. ACM Trans. Graph. 30, 6 (2011), 136.
  • [ZB05] Zhu Y., Bridson R.: Animating Sand as a Fluid. ACM Trans. Graph. 24, 3 (2005), 965–972.