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 largescale behavior. It is almost impossible for artists to add smallscale details to a lowresolution simulation without changing the largescale motion. Even more challenging is to modify the domain or rerun a simulation with different fluid parameters without affecting the largescale 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 divergencefree, 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 solidwall 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 nonseparating cells. We then enforce solidwall BCs for the velocity at nonseparating 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 firstorder PrimalDual (PD) method proposed by Pock et al. [PCBC09]. It features an optimal convergence rate for nonsmooth convex problems. Note that a whole class of primaldual 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 gridparticle 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 timeevolved velocity is divergencefree (i.e., masspreserving), which is typically implemented via Chorinlike pressureprojections, 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 positionbased 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.
Guiding
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 finetune, 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 highresolution simulations distant to the interface but are limited in volumetric contexts. Other approaches aiming for highlevel 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 lowfrequencies and pressureprojection. 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 multiscale approach based on control particles that controls lowspatial frequencies. Possibly the most similar approach to ours is the one by Nielsen et al. [NCZ09, NC10] who use a multiscale 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 nonsmooth optimization. Furthermore, for practical purposes, we introduce an upres method into our workflow [KTJG08, YCZ11, HMK11, RLL13, HK13] to facilitate the separation of lowresolution and highresolution 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 FFTbased one of R. Henderson [Hen12] are similar in spirit to our method, as they separate the BCs from the divergencefree 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
(1)  
(2) 
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 divergencefree condition in Eq. (2). A pressure field that exactly counteracts the divergence is computed to correct the velocity field. Orthogonality of the curlfree and divergencefree components ensures that divergencefree components of the flow are not affected and allows the pressure projection to be interpreted as an Euclidean projection onto the space of divergencefree velocity fields.
The splitting approach closely resembles convex optimization approaches originally developed for imaging inverse problems and machine learning involving nonsmooth 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
(3) 
where is a convex function that may be nonsmooth or even discontinuous. If is the sum of two simpler functions (say ), then we have
(4) 
A number of recently developed algorithms target this type of problem by employing an iterative divideandconquer approach. These algorithms are known as proximal methods and are defined in terms of socalled proximal operators (we use exclusively to denote the generic argument variable):
(5) 
One such algorithm is the PD method [PCBC09], which solves a slightly more general problem:
(6) 
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
(7)  
(8)  
(9) 
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:
(10) 
The variable updates are thus reduced to
(11)  
(12)  
(13) 
A more indepth 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 PDbased 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
(14)  
subject to 
where is the velocity field we seek, and is the space of divergencefree 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 divergencefree. 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
(15) 
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 periteration 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.
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 setup 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 periteration 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 divergencefree 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 largescale motions of a given target velocity. The objective function is given by
(16) 
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
(17) 
where
(18)  
(19)  
(20) 
Note that since is symmetric.
We can combine (three equations per cell in 3D) and the divergencefree constraint (one equation per cell) into a linear system . But since this system is overconstrained, we solve it in the leastsquares sense by considering
(21) 
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
(22) 
Factoring out , we obtain
(23) 
where
(24)  
(25) 
Computing for every iteration is slow. Instead, we precompute and since they do not rely on previous iterations. In fact, using the ShermanMorrisonWoodbury formula, can be approximated as
(26) 
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 straightforward 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.
Evaluation
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 . Detailpreserving guiding [TKRP06] subtracts the largescale (i.e., blurred) motions from the current velocity before adding the target velocity: . The idea is to create a new velocity with largescale motions from and smallscale 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)).
Detailpreserving blend allows guiding of largescale 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 largescale guiding strength and the blur radius controlling the smallscale 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.
Performance
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 fluidspecific 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 divergencefree neighbor but does not yield the divergencefree 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 E51650 CPUs.
For fairness of comparison, we experimented with several test scenes to deduce asoptimalaspossible 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.
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 
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 PrimalDual 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 noniterative method to greatly speed up the simulation. The method involves defining in terms of the nondivergent components of the input velocities and then finding its minimizer. Since differentiation commutes with convolution, the result will be automatically divergencefree for spatially invariant weights. Although this noniterative method offers potential speedup, we focus on the general case of spatially varying operators that requires iterating.
4.1 Guiding Results
To realize interesting guided flows, we found an iterative upres workflow that works well in practice. This is in line with previous work, where a lowresolution 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 starshaped velocity field with and . Then we run highresolution simulations with various values. If the same guiding is done without upres, 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 upres process for 3D simulations. Figure 10 shows a simple plume example. We first run a lowresolution () 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 smallscale 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 starshaped 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 lowresolution () 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.
Finally, Figure 11 demonstrates an example with obstacle, for which we upsampled from to with and . This highresolution version captures the input motion, but develops many interesting smallscale 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 lowresolution 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 Fourierdomain 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.
5 Separating SolidWall BCs
Another application of our PDbased method is the realization of separating solidwall 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 largescale fluid flow, the separation is preferable due to its improved realism.
We achieve this behavior by implementing separating solidwall 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 fluidsolid 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.
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 nonseparating 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: nonseparating (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 nonseparating, 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.
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 nonseparating cells based on the velocity field . The normal velocity components at nonseparating 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).
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 solidwall 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 fluidsolid 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 nonseparating 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 nonseparating 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 pseudocode 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 highresolution simulations.^{1}^{1}1 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 nonseparating 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 nondivergence 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 PDbased method performs on par with IOP, which we attribute to the low number of iterations required in this case; the higherorder convergence of the PD method does not develop its full potential in such a setting.
Figure 17 highlights that regular nonseparating 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.
6 Discussion and Conclusion
We presented a general framework for incorporating the PrimalDual method into fluid simulation, and demonstrated two applications: fluid guiding and separating BCs. We proposed a generic version of the PrimalDual based optimization scheme with fast convergence for general proximal operator subproblems. In addition, we discussed several extensions that are particularly wellsuited 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.
Limitations
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.
Outlook
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 highlevel flow guiding, into existing solvers without the need for complicated extensions.
Furthermore, there is a large number of interesting avenues to be explored with highlevel 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 (StG2015637014) 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 (applicationdependent)

: 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 nonseparating fluidsolid 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
(27)  
(28) 
The Krylov method can improve its convergence rate:
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 offdiagonal terms. By the ShermanMorrisonWoodbury Formula,
(32) 
Now since and both contains small value entries, is approximately zero. Hence
(33)  
(34) 
Note that the calculation of is trivial since is diagonal.
Appendix D Accelerated BC Solver
Below, we summarize our accelerated solver for separating solidwall BCs. Instead of fully handling the solidwall 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 solidwall 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 nonseparating cells (Neumann BCs for the pressure solver). We then iteratively set the velocity BCs, apply a regular CG pressure solver and reclassify 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 nonseparating 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 "".
References
 [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 solidfluid coupling. ACM Trans. Graph. 26, 3 (July 2007).
 [BHN07] Bridson R., Houriham J., Nordenstam M.: Curlnoise 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 NavierStokes equations. Mathematics of computation 22, 104 (1968), 745–762.
 [CMF12] Chentanez N., MuellerFischer 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 firstorder primaldual 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.: Targetdriven smoke animation. In ACM Transactions on Graphics (TOG) (2004), vol. 23, ACM, pp. 441–448.
 [GB13] Gerszewski D., Bargteil A. W.: Physicsbased animation of largescale 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, 34 (2000), 393–415.
 [HDN16] Heide F., Diamond S., Nießner M., RaganKelley 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.: Previewbased 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.: Highquality 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.: Freeflowing granular materials with twoway 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.: Primaldual 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 mumfordshah 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.: DetailPreserving Fluid Control. In Proceedings of the Symposium on Computer animation (2006), pp. 7–12.
 [YCZ11] Yuan Z., Chen F., Zhao Y.: Patternguided 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.