A Hybrid Lagrangian/Eulerian Collocated Advection and Projection Method for Fluid Simulation
Steven W. Gagniere David A. B. Hyde Alan Marquez-Razon Chenfanfu Jiang Ziheng Ge Xuchen Han Qi Guo Joseph Teran
cs.GR cs.CE physics.comp-ph








  UCLA
sgagniere@math.ucla.edu



  UCLA
dabh@math.ucla.edu



  UCLA
marqueza04@g.ucla.edu



  University of Pennsylvania
cffjiang@seas.upenn.edu



  UCLA
zihengge@math.ucla.edu



  UCLA
xhan0619@ucla.edu



  UCLA
qiguo@ucla.edu



  UCLA
jteran@math.ucla.edu





    < g r a p h i c s >

We simulate detailed incompressible flows with free surfaces and irregular domains using our novel hybrid particle/grid simulation approach.  Our numerical method yields intricate flow details with little dissipation, even at modest spatial resolution. Furthermore, we use collocated rather than staggered MAC grids.



We present a hybrid particle/grid approach for simulating incompressible fluids on collocated velocity grids. We interchangeably use particle and grid representations of transported quantities to balance efficiency and accuracy. A novel Backward Semi-Lagrangian method is derived to improve accuracy of grid based advection. Our approach utilizes the implicit formula associated with solutions of Burgers' equation. We solve this equation using Newton's method enabled by C^1 continuous grid interpolation. We enforce incompressibility over collocated, rather than staggered grids. Our projection technique is variational and designed for B-spline interpolation over regular grids where multiquadratic interpolation is used for velocity and multilinear interpolation for pressure. Despite our use of regular grids, we extend the variational technique to allow for cut-cell definition of irregular flow domains for both Dirichlet and free surface boundary conditions.









<ccs2012>
<concept>
<concept_id>10002950.10003714.10003715.10003750</concept_id>
<concept_desc>Mathematics of computing Discretization</concept_desc>
<concept_significance>500</concept_significance>
</concept>
<concept>
<concept_id>10002950.10003714.10003727.10003729</concept_id>
<concept_desc>Mathematics of computing Partial differential equations</concept_desc>
<concept_significance>500</concept_significance>
</concept>
<concept>
<concept_id>10002950.10003705.10003707</concept_id>
<concept_desc>Mathematics of computing Solvers</concept_desc>
<concept_significance>300</concept_significance>
</concept>
<concept>
<concept_id>10010405.10010432.10010441</concept_id>
<concept_desc>Applied computing Physics</concept_desc>
<concept_significance>300</concept_significance>
</concept>
</ccs2012>


[500]Mathematics of computing Discretization
[500]Mathematics of computing Partial differential equations
[300]Mathematics of computing Solvers
[300]Applied computing Physics










A Hybrid Lagrangian/Eulerian Collocated Advection and Projection Method for Fluid Simulation
    Joseph Teran
    December 30, 2023
============================================================================================





§ INTRODUCTION




Whether it be billowing smoke, energetic explosions, or breaking waves, simulation of incompressible flow is an indispensable tool for modern visual effects. Ever since the pioneering works of Foster and Metaxas foster:1996:liquids, Stam stam:1999:stable and Fedkiw et al. fedkiw:2001:visual,foster:2001:practical, the Chorin chorin:1967:numerical splitting of advective and pressure projection terms  has been the standard in computer graphics applications <cit.>. Most techniques use regular grids of Marker-And-Cell (MAC) <cit.> type with pressure and velocity components staggered at cell centers and faces respectively. Furthermore, advection is most often discretized using semi-Lagrangian techniques originally developed in the atmospheric sciences <cit.>. Although well-established, these techniques are not without their drawbacks. For example, the staggering utilized in the MAC grids is cumbersome since variables effectively live on four different grids. This can complicate many algorithms related to incompressible flow. E.g. Particle-In-Cell (PIC) <cit.> techniques like FLIP <cit.>, Affine/Polynomial Particle-In-Cell (APIC/PolyPIC) <cit.> and the Material Point Method (MPM) <cit.> must transfer information separately to and from each individual grid. Similarly, semi-Lagrangian techniques must separately solve for upwind locations at points on each of the velocity component grids. Moreover, while semi-Lagrangian techniques are renowned for the large time steps they admit (notably larger than the Courant-Friedrichs-Levy (CFL) condition), their inherent stability is plagued by dissipation that must be removed for most visual effects phenomena. Another limitation of the MAC grid arises with free-surface water simulation. In this case, the staggering prevents many velocity components near the fluid free surface from receiving a correction during projection (see e.g. <cit.>). Each of these velocity components must then be separately extrapolated to from the interior to receive a pressure correction. 



MAC grids are useful because the staggering prevents pressure null modes while allowing for accurate second order central differencing in discrete grad/div operators. However, there are alternatives in the computational physics literature. Many mixed Finite Element Method (FEM) techniques use collocated velocities <cit.> without suffering from pressure mode instabilities. For example, Taylor-Hood elements <cit.> use collocated multi-quadratic velocity interpolation and multilinear pressure interpolation to enforce incompressiblity. Recently, B-spline interpolation <cit.> has been used with Taylor-Hood <cit.>. We build on this work and develop an approach based on collocated multi-quadratic B-spline interpolation for velocities. This choice is motivated by the simplicity of collocated grids compared to staggering, but also because of the ease of attaining continuous derivatives with B-spline interpolation. For example, this interpolation is often chosen with MPM applications since C^1 interpolation is essential for stability <cit.>. In the context of fluids, we show that this allows for extremely stable and accurate advection.




We develop a new approach for Chorin splitting chorin:1967:numerical based on the collocated multiquadratic B-spline velocity, multilinear pressure Taylor-Hood element <cit.>. However, unlike the fully collocated technique of Bressan bressan:2010:isogeometric, we stagger pressures on the nodes of the grid and velocities at cell centers as in <cit.>, since it reduces coupling in the pressure projection system and naturally accommodates particle-based definition of the flow domain for free-surface simulation of water. Notably, our formulation does not require velocity extrapolation after pressure projection for free-surface flow calculations as is typically needed with MAC grids. We use regular grids, but as in <cit.>, we allow for irregular domains in a variational way using cut cells. However, rather than a weighted finite different approach, we use an FEM approach as in XFEM <cit.> and virtual node (VNA) <cit.> techniques. In VNA and XFEM approaches, integrals arising in the variational formulation are carried out over the intersection of the grid with the domain geometry. 




We leverage C^1 continuity guaranteed by our quadratic B-spline velocity interpolation to develop BSLQB, a novel Backward Semi-Lagrangian (BSL) <cit.> technique that achieves second order accuracy in space and time. BSL techniques utilize the implicit form of semi-Lagrangian advection. We show that our novel BSL method for quadratic B-splines dramatically reduces numerical dissipation with only a small modification to the widely-adopted explicit semi-Lagrangian formulations typically used in graphics applications. Semi-Lagrangian techniques for velocity advection utilize the implicit relation associated with solution of Burgers' equation

    (,t)=(-(t-s)(,t),s) D/Dt=∂/∂ t+∂/∂=0

for s≤ t <cit.>. Traditionally, graphics applications have preferred the explicit variant of semi-Lagragian advection whereby grid velocities are updated through the expression

    _^n+1=(_-Δ t _^n,t^n)

where _ is the location of grid node , _^n,_^n+1 are velocities at the node at times t^n and t^n+1 respectively and interpolation over the velocity grid is used to estimate (_-Δ t _^n,t^n) at non-grid node locations <cit.>. In contrast, BSL techniques leverage Equation  (<ref>) directly 

    _^n+1=(_-Δ t _^n+1,t^n)

which requires the solution of an implicit equation for _^n+1 <cit.>. Since our grid interpolation is naturally C^1, we show that this can be done very efficiently using a few steps of Newton's method. While this is more expensive than the explicit semi-Lagrangian formulations, we note that each node can still be updated in parallel since the implicit equations for _^n+1 are decoupled in . We show that solution of the implicit Equation (<ref>), rather than the traditionally used explicit Equation (<ref>) improves the order of convergence from first to second (in space and time). Notably, this does not require use of multiple time steps for backward/forward estimations of error, as is commonly done <cit.>. Furthermore, our method allows for larger-than-CFL time steps and is as stable or more so than explicit semi-Lagrangian formulations.




Lastly, we develop a hybrid particle/BSLQB advection technique that utilizes PolyPIC <cit.> in portions of the domain covered by particles and BSLQB in portions without particles. Our formulation naturally leverages the strengths of both approaches. Dense concentrations of particles can be added to regions of the domain where more detail is desired. Also, if particle coverage becomes too sparse because of turbulent flows, BSLQB can be used in the gaps. We demonstrate the efficacy of this technique with smoke simulation and narrow banding of particles near the fluid surface with water simulations as in <cit.>. In this case, level set advection naturally enabled with our BSLQB formulation is preferred in deeper water regions. We summarize our contributions as:


  * A novel cut-cell collocated velocity B-spline mixed FEM method for Chorin chorin:1967:numerical splitting discretization of the incompressible Euler equations.

  * BSLQB: a novel BSL technique designed for collocated multiquadratic B-spline velocity interpolation that achieves second order accuracy in space and time.

  * A hybrid BSLQB/PolyPIC method for narrow band free-surface flow simulations and concentrated-detail smoke simulations.






§ PREVIOUS WORK




 §.§ Advection

Stam stam:1999:stable first demonstrated the efficacy of semi-Lagrangian techniques for graphics applications and they have since become the standard, largely due to the large time steps they engender and their simple interpolatory nature. Many modifications to the original approach of Stam stam:1999:stable have been developed, often inspired by approaches in the engineering literature. Fedkiw et al. fedkiw:2001:visual use vorticity confinement <cit.> to counterbalance vorticity lost to dissipation and cubic grid interpolation. Kim et al. kim:2006:advections,kim:2005:bfecc and Selle et al. <cit.> combine forward and backward semi-Lagrangian steps to estimate and remove dissipative errors. Constrained Interpolation Profile <cit.> techniques additionally advect function derivatives to reduce dissipation. Molemaker et al. molemaker:2008:low use the QUICK technique of Leonhard leonard:1979:stable which is essentially upwinding with quadratic interpolation and Adams-Bashforth temporal discretization, although this does not have the favorable stability properties of semi-Lagrangian. Backward Difference Formula techniques are useful because they use an implicit multistep formulation for higher-order semi-Lagrangian advection yet still only require one projection per time step <cit.>.



The main idea in semi-Lagrangian techniques is to interpolate data from a characteristic point. This idea goes back to the Courant-Issaacson-Rees courant:1952:solution method. However, as noted in <cit.> semi-Lagrangian advection is very popular in atmospheric science simulation and the variants used in graphics that account for characteristics traveling beyond the local cell in one time step go back to Sawyer sawyer:1963:semi. The first BSL approach utilizing Equation (<ref>) was done by Robert robert:1981:stable in which they use fixed point iteration to solve the nonlinear equation. They fit a bicubic function to their data over 4×4 grid patches, then use that function in the fixed point iteration. If the upwind point leaves the grid, they clamp it to the boundary of the 4×4 patch. This clamping will degrade accuracy for larger time steps. In this case, more general interpolation is typically used (see <cit.> for useful reviews). Pudykiewicz and Staniforth pudykiewicz:1984:some investigate the effects of BSL versus explicit semi-Lagrangian. Specifically, they compare Bates and McDonald bates:1982:multiply (explicit) versus Robert robert:1981:stable (BSL). They show that keeping all things equal, the choice of Equation (<ref>) (explicit) instead of Equation (<ref>) (BSL) leads to more dissipation and mass loss. This is consistent with our observations with BSLQB.



Interestingly, multiquadratic B-splines have not been adopted by the semi-Lagrangian community, despite their natural regularity. Hermite splines, multicubic splines and even Lagrange polynomials are commonly used <cit.>. Preference for Hermite splines and Lagrange polynomials is likely due to their local nature (they do not require solution of a global system for coefficients) and preference for multicubic splines (over multi-quadratic) is possibly due to the requirement of odd degree for natural splines (odd degree splines behave like low pass filters and tend to be smoother than even degree splines <cit.>). Cubic splines are considered to be more accurate than Hernite splines and Lagrange interpolation <cit.>. Interestingly, Riishøjgaard et al. riishojgaard:1998:use found that cubic spline interpolation gave rise to a noisier solution than cubic Lagrange interpolation with a technique analogous to that of Makar and Karpik makar:1996:basis. However, they also note that addition of a selective scale diffusion term helps reduce noise associated with cubic splines. Wang and Layton wang:2010:new use linear B-splines with BSL but only consider one space dimension which makes Equation (<ref>) linear and easily solvable.



Dissipation with explicit semi-Lagrangian advection is so severe that many graphics researchers have resorted to alternative methods to avoid it. Mullen et al. mullen:2009:energy develop energy preserving integration to prevent the need for correcting dissipative behavior. Some authors <cit.> resolve the flow map characteristics for periods longer than a single time step (as opposed to one step with semi-Lagrangian) to reduce dissipation. Hybrid Lagrange/Eulerian techniques like PIC (and related approaches) <cit.> explicitly track motion of particles in the fluid, which is nearly dissipation-free, but can suffer from distortion in particle sampling quality. Vorticity formulations are also typically less dissipative, but can have issues with boundary conditions enforcement <cit.>. Zehnder et al., Zhang et al. and Mullen et al. mullen:2009:energy,zehnder:2018:advection,narain:2019:ref,zhang:2015:restoring have noted that the Chorin projection itself causes dissipation. Zhang et al. zhang:2015:restoring reduced artificial dissipation caused by the projection step by estimating lost vorticity and adding it back into the fluid. Zehnder et al. zehnder:2018:advection,narain:2019:ref propose a simple, but very effective modification to the splitting scheme that is similar to midpoint rule integration to reduce the projection error.



 §.§ Pressure projection


Graphics techniques utilizing pressure projection typically use voxelized MAC grids with boundary conditions enforced at cell centers and faces, however many methods improve this by taking into account sub-cell geometric detail. Enright et al. enright:2003:using showed that enforcing the pressure free surface boundary condition at MAC grid edge crossings (rather than at cell centers) dramatically improved the look of water surface waves and ripples. Batty, Bridson and colleagues developed variational weighted finite difference approaches to enforce velocity boundary conditions with MAC grids on edge crossings and improved pressure boundary conditions at the free surface in the case of viscous stress <cit.>. XFEM <cit.> and virtual node (VNA) <cit.> techniques also use cut cell geometry with variational techniques. Schroeder et al. schroeder:2014:vna use cut cells with MAC grids, but their technique is limited to moderate Reynolds numbers. 




There is a vast literature on enforcing incompressibility in the FEM community <cit.>. Our approach is most similar to the B-spline Taylor-Hood element of Bressan <cit.>. Adoption of B-spline interpolation in FEM is part of the isogeometric movement <cit.>. Originally motivated by the desire to streamline the transition from computer-aided design (CAD) to FEM simulation, isogeometric analysis explores the use of CAD-based interpolation (e.g. B-splines and nonuniform rational B-splines (NURBS)) with FEM methodologies. Hughes et al. hughes:2005:isogeometric show that in addition to simplifying the transition from CAD to simulation, the higher regularity and spectral-like properties exhibited by these splines makes them more accurate than traditionally used interpolation. We enforce Dirichlet boundary conditions weakly as in XFEM and VNA approaches <cit.>. Bazilevs et al. bazilevs:2007:weak show that weak Dirichlet enforcement with isogeometric analysis can be more accurate than strong enforcement.



Graphics applications are typically concerned with turbulent, high-Reynolds numbers flows. Interestingly, B-splines have proven effective for these flows by researchers in the Large Eddy Simulation (LES) community <cit.>. Kravchenko et al. kravchenko:1999:bspline use a variational weighted residuals approach with B-splines for turbulent LES and show that the increased regularity significantly reduces computational costs. Boatela et al. botella:2002:collocation use a similar approach, but apply a collocation technique where the strong form of the div-grad formulation of incompressibility is enforced point wise. They show that their B-spline approach attains optimal order of accuracy with accurate resolution of quadratic flow invariants. Boatela et al. botella:2002:collocation also introduce a notion of sparse approximation to the inverse mass matrix to avoid dense systems of equations in the pressure solve.





§ GOVERNING EQUATIONS AND OPERATOR SPLITTING

We solve the incompressible Euler equations that describe the evolution of a fluid in terms of its mass density ρ, velocity , pressure p and gravitational constant ≫ as

    ρD/Dt   =ρ(∂/∂ t + ∂/∂)=-∇ p + ρ≫,  ∈Ω
    ∇·   = 0,  ∈Ω
    ·   =a,  ∈∂Ω_D 
    
    p   =0,  ∈∂Ω_N

where Equation (<ref>) is balance of linear momentum, Equation (<ref>) is the incompressibility constraint, Equation (<ref>) is the boundary condition for the normal component of the velocity and Equation (<ref>) is the free surface boundary condition. We use Ω to denote the region occupied by the fluid, ∂Ω_D to denote the portion of the boundary of the fluid domain on which velocity is prescribed to be a (which may vary over the boundary) and ∂Ω_N is the surface of the water where the pressure is zero (see Figure <ref>).




In a Chorin chorin:1967:numerical operator splitting of the advective and pressure terms, velocity is first updated to an intermediate field  under the convective ρD/Dt=0, followed by an update from the pressure and gravitational body forcing under ρ∂/∂ t=-∇ p + ρ≫ where the pressure is determined to enforce ∇· = 0. Dividing by the mass density, the convective step is seen to be an update under Burgers' equation (<ref>). Burgers' equation governs temporally constant Lagrangian velocity (zero Lagrangian acceleration). The characteristic curves for flows of this type are straight lines (since the Lagrangian acceleration is zero), on which the velocity is constant (see Figure <ref>). This gives rise to the implicit relation (,t)=(-(t-s)(,t),s) for s≤ t. Intuitively, if we want to know the velocity (,t) at point  at time t, we look back along the characteristic passing through  at time t to any previous time s; however, the characteristic is the straight line defined by the velocity (,t) that we want to know. Hence we take an implicit approach to the solution of this equation, which when combined with the operator splitting amounts to

    -^n/Δ t   =0
    ρ^n+1-/Δ t   =-∇ p^n+1 + ρ≫
    ∇·^n+1   = 0

where we use the notation ^n+α()=(,t^n+α), α=0,1 to denote the time t^n+α velocities. Furthermore, the intermediate velocity  is related to ^n through ^n()=(-Δ t (),t^n).





§ SPATIAL DISCRETIZATION

We discretize in space by first representing velocity and pressure in terms of mulitquadratic and multilinear B-splines for velocity and pressure respectively. We use a regular grid with spacing Δ x and define pressure degrees of freedom at grid vertices and velocity degrees of freedom at grid cell centers as in <cit.> (see Figure <ref>). This efficiently aligns the support of the multiquadratic and multilinear interpolating functions which naturally allows for a grid-cell-wise definition of the flow domain (see Figure <ref>). We use N_() to represent the multiquadratic B-spline basis function associated with velocity degree of freedom _ at grid cell center _ and χ_() for the multilinear basis function associated with pressure p_ at grid node _. These are defined as

    N_()   =∏_αN̂(x_α-x_α/Δ x),  χ_()=∏_αχ̂(x_α-x_α/Δ x)
    N̂(η)   ={[ (η +3/2)^2/2, η∈(-3/2,-1/2);     -η^2+3/4,  η∈[-1/2,1/2]; (η -3/2)^2/2,   η∈(1/2,3/2);            0,     otherwise ].  χ̂(ν)={[      1+ν,  ν∈(-1,0);      1-ν,   ν∈[0,1);        0, otherwise ].

where we use Greek indices α to indicate components of the vectors , _ and _. With this convention we interpolate to define velocity and pressure fields 

    ()=∑__ N_(),  p()=∑_ p_χ_().

We use the notation _ to distinguish it from the velocity at the grid node (_)=∑__ N_(_) since the multiquadratic B-splines are not interpolatory and these will in general be different. Note that multilinear interpolation is interpolatory and p_=∑_ p_χ_(_).


 §.§ BSLQB Advection

With this interpolation choice, we first solve for intermediate grid node velocity values (_) from Equation (<ref>) as

    (_)=∑_^n_ N_(_-Δ t (_)).

We can solve this equation using Newton's method since the multiquadratic B-splines are C^1. We use ^k_ to denote the k^th Newton approximation to (_). Explicit semi-Lagrangian is used as an initial guess with ^0_=∑_^n_ N_(_-Δ t ∑_≪^n_≪ N_≪(_)) and then we update iteratively via ^k_+=δ^k with Newton increment δ^k satisfying

    δ^k   =(+Δ t ∂^n/∂(_-Δ t ^k_))^-1(∑_^n_ N_(_-Δ t ^k_)-^k_)

where ∂^n/∂(_-Δ t ^k_)=∑_^n_∂ N_/∂(_-Δ t ^k_). It is generally observed <cit.> that with BSL approaches of this type, this iteration will converge as long as +Δ t ∑_^n_∂ N_/∂(_-Δ t ^k_) is non-singular. We note that this condition holds as long as no shocks form under Burgers' equation <cit.> (forward from time t^n). This is a safe assumption since we are modeling incompressible flow with which shock formation does not occur, but it may be a problem for compressible flows. In practice, this iteration converges in 3 or 4 iterations, even with CFL numbers larger than 4 (see Section <ref>). When it does fail (which occurs less than one percent of the time in the examples we run), it is usually for points near the boundary with characteristics that leave the domain (since we cannot estimate ∂^n/∂ using grid interpolation if the upwind estimate leaves the grid). In this case we use explicit semi-Lagrangian and interpolate from the boundary conditions if the characteristic point is off the domain.



Once we have obtained the grid node values of the intermediate velocity (_), we must determine interpolation coefficients _ such that (_)=∑__ N_(_). On the boundary of the grid, we set _ = (_) since we can only interpolate to _ if all of its neighbors have data. This yields a square, symmetric positive definite system of equations for the remaining _. The system is very well conditioned with sparse, symmetric matrix N_(_) consisting of non-negative entries and rows that sum to one.  The sparsity and symmetry of the system arises from the compact support and geometric symmetry, respectively, of the B-spline basis functions N_. The system can be solved to a residual of machine precision in one iteration of PCG (or tens of iterations of unpreconditioned CG). In practice, we have noticed that for some flows, determining the coefficients _ can lead to increasingly oscillatory velocity fields. This is perhaps due to the unfavorable filtering properties of even order B-splines  <cit.>. However, we found that a simple stabilization strategy can be obtained as

    ∑_(λ N_(_) + (1-λ)δ_)_=(_)

where λ∈[0,1] and δ_ is the Kronecker delta. A value of λ=0 is very stable, but extremely dissipative. Stable yet energetic behavior is achieved by decreasing the value of λ under grid refinement. In practice we found that λ∈ (.95,1 ] with λ =cΔ x for constant c provided a good balance without compromising second order accuracy of the method (see Section <ref>). We note that Riishøjgaard et al. riishojgaard:1998:use also added diffusion to cubic spline interpolation based semi-Lagrangian to reduce noise.


 §.§ Hybrid BSLQB-PolyPIC Advection

In some portions of the domain, we store particles with positions _p^n and PolyPIC <cit.> velocity coefficients ^n_p. In the vicinity of the particles, we use PolyPIC <cit.> to update the intermediate velocity field _. First we update particle positions as ^n+1_p=^n_p + Δ t _p^n (where the velocity _p^n is determined from ^n_p following <cit.>). Then the components w̅_α of the coefficients _ are determined as

    w̅_α=∑_p m_pN_(_p^n+1)(∑_r=1^N_r s_r(_-^n+1_p)c^n_prα)/∑_p m_pN_(_p^n+1)

where N_r is the number of polynomial modes s_r(), as in Fu et al. fu:2017:poly. To create our hybrid approach, we update  w̅_α from Equation (<ref>) whenever the denominator is greater than a threshold ∑_p m_pN_(_p^n+1)>τ^m, otherwise we use the BSLQB update from Equation (<ref>). We use this threshold because the grid node update in Equation(<ref>) loses accuracy when the denominator is near zero and in this case the BSLQB approximation is likely more accurate. Note that the polynomial mode coefficients for the next time step ^n+1_p are determined from the grid velocities at the end of the time step (using particle positions _p^n+1 and after pressure projection).



§ PRESSURE PROJECTION

We solve Equations (<ref>)-(<ref>) and boundary condition Equations (<ref>)-(<ref>) in a variational way. To do this, we require that the dot products of Equations (<ref>), (<ref>) and Equations (<ref>) with arbitrary test functions , q and μ respectively integrated over the domain are always equal to zero. The free surface boundary condition in Equation (<ref>) is naturally satisfied by our treatment of Equation (<ref>). We summarize this as 

    ∫_Ω·ρ(^n+1-/Δ t)d   =∫_Ω p^n+1∇· + ρ·≫ d
       -∫_∂Ωp^n+1· ds()
    ∫_Ω q ∇·^n+1 d   =0 
    ∫_∂Ω_Dμ(^n+1·-a) ds()   =0.

Here we integrate by parts in the integral associated with Equation (<ref>). Furthermore, we modify the expression ∫_∂Ωp^n+1· ds() in Equation (<ref>) in accordance with the boundary conditions. We know that the pressure is zero on ∂Ω_N, however we do not know its value on ∂Ω_D. We introduce the pressure on this portion of the domain as a Lagrange multiplier λ^n+1 associated with satisfaction of the velocity boundary condition in Equation (<ref>). Physically, this is the external pressure we would need to apply on ∂Ω_D to ensure that ^n+1·=a. With this convention, we have ∫_∂Ωp^n+1· ds()=∫_∂Ω_Dλ^n+1· ds(). We note that unlike Equation (<ref>) (and its strong form counterpart (<ref>)) that requires introduction of a Lagrange multiplier, Equation (<ref>) is naturally enforced through the weak form simply by setting p^n+1=0 in the integral over ∂Ω_N in Equation (<ref>).



To discretize in space, we introduce interpolation for the test functions , q and μ. We use the same spaces as in Equation (<ref>) for velocity and pressure for =∑__ N_ and q=∑_ q_χ_. For the test functions μ, we choose the same space as q,p, but with functions restricted to ∂Ω_D, μ=∑_μ_χ_ for  with grid cell Ω_∩∂Ω_D ≠∅ (see Figure <ref>). We choose the same space for λ^n+1=∑_λ^n+1_χ_ to close the system. With these choices for the test functions, the variational problem is projected to a finite dimensional problem defined by the interpolation degrees of freedom. This is expressed as a linear system for velocities _^n+1, internal pressures p^n+1_, and external pressures λ_^n+1 that is equivalent to 

    ([     -^T  ^T;   -;         ])
    (
    [  ^n+1;  ^n+1; Λ^n+1 ])=
    (
    [ + ≫̂;    0;    Å ]).

Here ^n+1, ^n+1 and Λ^n+1 are the vectors of all unknown _^n+1, p_^n+1 and λ_^n+1 respectively. Furthermore  is the mass matrix,  defines the velocity boundary conditions and  defines the discrete divergence condition. Lastly,  is the vector of all _ that define the intermediate velocity, ≫̂ is from gravity and Å is the variational boundary condition. Using the convention that Greek indices α,β range from 1-3, these matrices and vectors have entries 

    M_αβ=δ_αβ∫_Ωρ/Δ t N_ N_ d,  D_β   =∫_Ωχ_∂ N_/∂ x_βd,  ĝ_α=∫_Ωρ g_α N_ d
    
    B_β=∫_Ω_Dχ_ N_ n_β ds(),  A_   = ∫_Ω aχ_ ds().

If we define =[-^T,^T], we can convert this system into a symmetric positive definite one for ^n+1 and Λ^n+1 followed by a velocity correction for ^n+1

    (
    [  ^n+1; Λ^n+1 ])   =(^T^-1)^-1(^T(+^-1≫̂)-([ 0; Å ]))
    ^n+1   =-^-1(
    [  ^n+1; Λ^n+1 ])++^-1≫̂.

Unfortunately, this system will be dense in the current formulation since the full mass matrix M_αβ is non-diagonal with dense inverse <cit.>. However, a simple lumped mass approximation

    M^l_αβ={[ δ_αβ∫_Ωρ/Δ t N_ d,                  =;                 0,          otherwise ].

gives rise to a sparse matrix in Equation (<ref>).
























 §.§ Cut cells


As in XFEM and VNA approaches <cit.>, we resolve sub-grid-cell geometry by simply performing the integrations in Equations (<ref>)-(<ref>) over the geometry of the fluid domain. We use a level set to define solid boundaries (green in Figure <ref>) on which velocity boundary conditions are defined. We triangulate the zero isocontour using marching cubes <cit.> (see Figure <ref>). The integrals in Equations (<ref>)-(<ref>) all involve polynomials over volumetric polyhedra (Equations (<ref>), blue in Figure <ref>) or surface polygons (Equations (<ref>), green in Figure <ref>) and we use Gauss quadrature of order adapted to compute the integrals with no error (see <cit.>). For free surface flows, we use particles (and additionally a level set function in the case of narrow banding, see Section (<ref>)) to denote grid cells with fluid in them. Cells near the solid boundary are clipped by the marching cubes geometry. The fluid domain Ω is defined as the union of all clipped and full fluid cells (see Figure <ref>). 
 


Notably, taking a cut cell approach with our variational formulation allows us to prove that our method can resolve a standing pool of water exactly without producing numerical currents. We know that with gravitational force ρ≫ (e.g. with ≫ pointing in the y direction with magnitude g), steady state is maintained if the pressure increases with depth as p=ρ g(y_0-y) where y_0 is the height of the water surface at rest, since -∇ p + ρ≫ = 0. Since we use multilinear interpolating functions for p, the exact solution is representable in our discrete space and with a short proof we show (see <cit.>) that this means our method will choose it to maintain a standing pool of water, independent of fluid domain boundary geometry.





§ NARROW BAND FREE SURFACE

For free surface flows, we develop a narrow band approach as in chentanez:2015:coupling,ferstl:2016:narrow,sato:2018:nb. We represent the fluid domain with a level set and seed particles in a band of width W from the zero isocontour (see Figure <ref>). Particles are advected and used to augment BSLQB advection as detailed in Section <ref>. We also advect the level set by interpolating its value at the previous step from the upwind location _-Δ t (_) determined in Equation (<ref>). We then use the updated particle locations to compute a narrow band level set from the particles based on the method of Boyd and Bridson <cit.>. We update the level set to be the union of that defined by the narrow band and that from advection. This is done by taking the minimum of the two level set values and then redistancing with the method of Zhao zhao:2005:fast.


§ EXAMPLES







 §.§ Hybrid BSLQB/PolyPIC

We demonstrate our hybrid BSLQB/PolyPIC advection with water simulation. We prevent excessive run times by utilizing a narrow band of particles near the free surface and a level set (with BSLQB advection) in deeper levels. Figure <ref> Top shows a disc of water splashing in a rectangular tank with dimension 1×2 and grid cell size Δ x = 1/255. The time step is restricted to be in the range Δ t ∈[0.005, 0.01].  20 particles are initialized in every cell that is initially in a narrow band of 7 Δ x below the zero isocontour of the level set. Figure <ref> Bottom shows an analogous 3D example where a sphere of water splashes in a tank.  A cell size of Δ x = 1/63 is used in a domain with dimensions 1× 2 × 1.  We take a fixed time step of Δ t = 0.01 and demonstrate that narrow banding does not prevent larger-than-CFL time steps.  1,008,187 particles are used to resolve the free surface in a narrow band of width 5Δ x.  As in 2D, the particles capture highly-dynamic behavior of the free surface while the level set is sufficient to represent the bulk fluid in the bottom half of the domain.



We also demonstrate our hybrid advection with a vortex shedding example (see Figure <ref>). The flow domain Ω is a 3× 1 rectangle with circle of radius 0.05. We seed a band of particles of width .2 above the midline y=.5 for PolyPIC advection. Advection in the rest of the domain is done with BSLQB. The vorticity plot illustrates a seamless transition between the two advection schemes. The simulation was run with a grid resolution of Δ x=1/255, CFL number of 4 (i.e. Δ t = 4Δ x/v_max), and inlet speed of 1.5.




 §.§ BSLQB comparison with explicit semi-Lagrangian

We demonstrate improved resolution of flow detail with BSLQB compared to explicit semi-Lagrangian in a 2D example of smoke flowing past a circle (see Figure <ref>) and with a 2D spinning circle example (see Figure <ref>). Note that particles are only used for flow visualization and not for PolyPIC advection in these examples. BSLQB exhibits more energetic, turbulent flows than semi-Lagrangian advection. Notably, the BSLQB result breaks symmetry sooner. In Figure <ref> we also examine the effect of extremal values of the λ parameter described in Equation (<ref>).  A zero value of λ is quite dissipative compared to a full value of λ = 1 for both semi-Lagrangian and BSLQB.  As mentioned in Section <ref>, we generally found that keeping λ close to 1 provided the least dissipative behavior, while setting the value slightly less than 1 helped restore stability when necessary (one can also dynamically adjust this value over the course of a simulation). In Figure <ref>, we initially set the angular velocity to 4 radians per second in a circle of radius .2 (with Ω=[0,1]×[0,1]). The simulation is run with Δ x=1/511 and a Δ t = .02 (CFL number of 3).



We examine the convergence behavior of BSLQB for the 2D Burgers' equation D/Dt=0 with initial data ()=·(Å) for Å=Λ^T for diagonal Λ with entries 1 and .25 and rotation (of .1 radians)  (see Figure <ref>). We examine the convergence behavior under refinement in space and time with Δ t=Δ x. We compute the best fit line to the plot of the logarithm of the L^∞ norm of the error versus the logarithm of Δ x for a number of grid resolutions. We observe slopes of approximately 2 for BSLQB with interpolation parameter λ=1 and λ=1-cΔ x (with c = 2.95), indicating second order accuracy in space and time under refinement. We observe slopes of approximately 1 for explicit semi-Lagrangian, indicating first order.







 §.§ Cut cell examples

We demonstrate the ability of our cut cell method to produce detailed flows in complicated irregular domains for smoke and free surface water examples. Figure <ref> demonstrates the subtle and visually interesting behavior that arises as two plumes of multicolored smoke flow to the center of a cubic domain colliding with a spherical boundary. We use Δ x = 1/63 and Δ t = .02. We demonstrate a more complex domain in Figure <ref>. Puffs of colored smoke with converging initial velocities are placed in a bunny shaped clear domain. We use grid size Δ x = 1/127 and a fixed time step of Δ t = 0.01 (CFL number >1). In Figure <ref>, we demonstrate water splashing, while accurately conforming to the walls of an irregular domain defined as the interior of a large sphere and exterior of a small inner sphere. The spatial resolution of the domain is Δ x = 1/127, and 30 particles per cell are seeded in the initial fluid shape. A minimum time step of Δ t=0.001 is enforced, which is often larger than the CFL condition. We also consider dam break simulations in rectangular domains with column obstacles (Figure <ref>) and a bunny obstacle (Figure <ref>). Both examples use a grid cell size of Δ x = 1/127, 8 particles per cell and a fixed time step of Δ t = 0.003. Lastly, we demonstrate the benefits of our cut cell formulation over a more simplified, voxelized approach in Figure <ref>. Notice the water naturally sliding in the cut cell domain compared with the jagged flow in the voxelized domain.


























 §.§ Performance considerations




The implementation of our method takes advantage of hybrid parallelism (MPI, OpenMP, and CUDA/OpenCL) on heterogeneous compute architectures in order to achieve practical runtime performance (see Table <ref> for 3D example performance numbers). The spatial domain is uniformly divided into subdomains assigned to distinct MPI ranks, which distributes much of the computational load at the expense of synchronization overhead exchanging ghost information across ranks. On each rank, steps of our time integration loop such as BSLQB advection are multithreaded using OpenMP or CUDA when appropriate. The dominant costs per time step are the solution of the pressure projection system and, in the case of free surface simulation, assembly of the pressure system and its preconditioner. We permute Equation (<ref>) so that each rank's degrees of freedom are contiguous in the solution vector then solve the system using AMGCL <cit.> using the multi-GPU VexCL backend (or the OpenMP CPU backend on more limited machines). Using a strong algebraic multigrid preconditioner with large-degree Chebyshev smoothing allows our system to be solved to desired tolerance in tens of iterations, even at fine spatial resolution.  An important step in minimizng the cost of system assembly is to scalably parallelize sparse matrix-matrix multiplication, for which we use the algorithm of Saad saad:2003:sparse.  In the future, we are interested in implementing load balancing strategies such as the simple speculative load balancing approach of <cit.>, particularly for free surface flows. We note that our implementation enables high-resolution simulations such as that in Figure <ref> at relatively modest computational cost (see Table <ref>).



§ DISCUSSION AND LIMITATIONS

Our approach has several key limitations that could be improved. First, our adoption of collocated multiquadratic velocity and multilinear pressure is a significant departure from most fluid solvers utilized in graphics applications. We note that BSLQB and BSLQB/PolyPIC could be used with a MAC grid; however, each velocity face component would have to be solved for individually. Another drawback for our multiquadratic velocity and multilinear pressure formulation is that it gives rise to a very wide pressure system stencil consisting of 49 non-zero entries per row in 2D and 343 in 3D. Collocated approaches that make use of multilinear velocities and constant pressure give rise to 9 (2D) and 27 (3D) entries per row <cit.>, however they do not allow for C^1 continuity and require spurious pressure mode damping. Our wide stencils likely negatively affect the efficacy of preconditioning techniques as well, however we were very pleased with the efficiency of the AMGCL <cit.> library. Also, while the use of mass lumping in Equation (<ref>) is necessary to ensure a sparse pressure projection system, Boatella et al. botella:2002:collocation note that this has been shown to degrade accuracy. In fact, Boatela et al. botella:2002:collocation introduce a sparse approximate inverse to the full mass matrix to avoid dense systems of equations in the pressure solve without degrading accuracy. Split cubic interpolation, which approximates similar systems with tridiagonal ones could also possibly be used for this <cit.>. Adoption of one of these approaches with our formulation would be an interesting area of future work. Also, we note that the more sophisticated transition criteria for narrow banding techniques in Sato et al. sato:2018:nb could naturally be used with our method.  Finally, we note that the work of Zehnder et al. zehnder:2018:advection,narain:2019:ref could be easily applied to our technique to further reduce dissipation since it is based on the Chorin chorin:1967:numerical splitting techniques (Equations (<ref>)-(<ref>)) that we start from.



ACM-Reference-Format


