@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Computer Science, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Zhang, Xinxin"@en ; dcterms:issued "2017-05-02T15:52:56Z"@en, "2017"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """Turbulent gaseous phenomena, often visually characterized by their swirling nature, are mostly dominated by the evolution of vorticity. Small scale vortex structures are essential to the look of smoke, fire, and related effects, whether produced by vortex interactions, jumps in density across interfaces (baroclinity), or viscous boundary layers. Classic Eulerian fluid solvers do not cost-effectively capture these small-scale features in large domains. Lagrangian vortex methods show great promise from turbulence modelling, but face significant challenges in handling boundary conditions, making them less attractive for computer graphics applications. This thesis proposes several novel solutions for the efficient simulation of small scale vortex features both in Eulerian and Lagrangian frameworks, extending robust Eulerian simulations with new algorithms inspired by Lagrangian vortex methods."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/61435?expand=metadata"@en ; skos:note "Tackling Small-Scale Fluid Features in Large Domainsfor Computer GraphicsbyXinxin ZhangB. Information and Computation Science, Zhejiang University City College, 2009M. Computer Science, New York University, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Computer Science)The University of British Columbia(Vancouver)April 2017c© Xinxin Zhang, 2017AbstractTurbulent gaseous phenomena, often visually characterized by their swirling na-ture, are mostly dominated by the evolution of vorticity. Small scale vortex struc-tures are essential to the look of smoke, fire, and related effects, whether pro-duced by vortex interactions, jumps in density across interfaces (baroclinity), orviscous boundary layers. Classic Eulerian fluid solvers do not cost-effectively cap-ture these small-scale features in large domains. Lagrangian vortex methods showgreat promise from turbulence modelling, but face significant challenges in han-dling boundary conditions, making them less attractive for computer graphics ap-plications. This thesis proposes several novel solutions for the efficient simulationof small scale vortex features both in Eulerian and Lagrangian frameworks, ex-tending robust Eulerian simulations with new algorithms inspired by Lagrangianvortex methods.iiPrefaceThis thesis incorporates three published papers as separate chapters:• Chapter 2: X. Zhang, M. Li, and R. Bridson. 2016. Resolving fluid boundarylayers with particle strength exchange and weak adaptivity. ACM Trans.Graph. 35,4(Proc. SIGGRAPH), Article 76.• Chapter 3: X. Zhang, R. Bridson, and C. Greif. 2015. Restoring the miss-ing vorticity in advection-projection fluid solvers, ACM Trans. Graph. 34, 4(Proc. SIGGRAPH), Article 52.• Chapter 4: X. Zhang and R. Bridson. 2014. A PPPM fast summation methodfor fluids and beyond. ACM Trans. Graph. 33, 6 (Proc. SIGGRAPH Asia),Article 206.Apart from the usual supervision roles by Dr. Robert Bridson and Dr. ChenGreif, the single-scattering particle renderer used for images and animations wasprovided by Dr. Bridson, and help from Minchen Li in testing simulation examplesfor Chapter 2, all work is mine.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Resolving Fluid Boundary Layers with Particle Strength Exchangeand Weak Adaptivity . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.3 VFLIP and weakly higher resolution regional projection (WHIRP). 102.3.1 Solving the convection-diffusion equation with ghost par-ticles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3.2 Regional projection for particle velocity correction . . . . 132.3.3 Our pressure solver . . . . . . . . . . . . . . . . . . . . . 142.3.4 Seeding and deleting particles . . . . . . . . . . . . . . . 152.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . 162.5 Conclusions, limitations, and future work . . . . . . . . . . . . . 20iv3 Restoring the Missing Vorticity in Advection Projection Fluid Solvers 233.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.2 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.3 The IVOCK scheme . . . . . . . . . . . . . . . . . . . . . . . . . 293.3.1 Vortex dynamics on grids . . . . . . . . . . . . . . . . . . 313.3.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 343.4 Applications and results . . . . . . . . . . . . . . . . . . . . . . . 373.4.1 Smoke . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.4.2 Liquids . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.4.3 Fire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424 A PPPM Fast Summation Method for Fluids and Beyond . . . . . . 444.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.2 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.3 Particle-particle particle-mesh method(PPPM) . . . . . . . . . . . 494.3.1 Particle-mesh step in open space . . . . . . . . . . . . . . 514.3.2 Cancelling local influences in the grid . . . . . . . . . . . 534.3.3 Velocity evaluation . . . . . . . . . . . . . . . . . . . . . 564.3.4 PPPM discussion . . . . . . . . . . . . . . . . . . . . . . 574.4 Vortex-Particle Smoke . . . . . . . . . . . . . . . . . . . . . . . 604.4.1 Vortex segment, vortex stretching and stability issue . . . 614.4.2 Solid boundaries and vortex shedding . . . . . . . . . . . 644.5 Results and conclusions . . . . . . . . . . . . . . . . . . . . . . . 674.5.1 PPPM for vortex flow . . . . . . . . . . . . . . . . . . . . 674.5.2 Applying PPPM to Poisson surface reconstruction . . . . 694.6 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.7 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 725 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77vList of TablesTable 2.1 Symbol abbreviations used throughout this paper. . . . . . . . 7Table 2.2 The Semi-AMGPCG solver applied to a Poisson problem withsingle impulse, in a unit cube domain with a zero Neumannboundary condition along a sphere with diameter 0.2 in the cen-ter, and zero Dirichlet boundary conditions at the edges of thecube. The number of iterations required to reduce the residualby a factor of 10−10 is essentially independent of the grid size,providing one extra digit of accuracy per iteration. . . . . . . . 15Table 2.3 Vortex shedding frequency statistics from numerical simulations 18Table 2.4 Timings of various simulations associated with this paper, mea-sured on a desktop machine with Intel(R) Core(TM) i7-3930KCPU and 32 GB RAM. . . . . . . . . . . . . . . . . . . . . . 20Table 2.5 Timings of flow past sphere simulations with only differencesin outer and inner grid sizes, measured on a desktop machinewith Intel(R) Core(TM) i7-3930K CPU and 32 GB RAM. . . . 20Table 3.1 Algorithm abbreviations used through out this paper. . . . . . . 23Table 3.2 Performance comparison of IVOCK augmenting different schemes,for a smoke simulation at 128x256x128 grid resolution, runningon an Intel(R) Core(TM) i7-3630QM 2.40GHz CPU. . . . . . 39Table 4.1 Common symbols used throughout the paper. . . . . . . . . . . 45Table 4.2 Accuracy of different method with and without the monopoleB.C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60viList of FiguresFigure 1.1 The chaotic motion of eddies produced by the storms in theatmosphere of Jupiter. Courtesy NASA/JPL-Caltech. . . . . . 2Figure 2.1 Simulation of a rotating fan in wind, Re ≈ 105, visualizedwith smoke marker particles. Large scale motion and detailedboundary layer dynamics are cheaply and easily coupled usingour method, producing high quality fluid animations at lowcost. Top: top view of fan wake at frames 80, 120 and 200.Bottom: side view at frame 300. . . . . . . . . . . . . . . . . 5Figure 2.2 Vortex shedding from a cylinder in 2D at Re=15000, zoom-in view near the boundary. Left column: simulation results atlow resolution (approximately 50×50 grid cells shown). Mid-dle column: simulation obtained with 4× resolution. Rightcolumn: simulation obtained with our method, with the samecoarse grid but a 4× refined grid overlaid just around the solid.Our method produces results visually consistent with the highresolution reference because only the boundary layer needsthat resolving power. . . . . . . . . . . . . . . . . . . . . . . 6Figure 2.3 Our method captures the velocity field near the boundary effi-ciently and with high apparent fidelity. Left: zoom-in near theboundary flow. Right: the entire simulation, Re=15000. . . . . 7viiFigure 2.4 A typical domain construction used in our method. Ω indicatesthe global Eulerian domain where fluid motion is loosely cap-tured. Eulerian subdomains with finer resolution can be placedanywhere to enhance the simulation quality locally, such as thegreen domain(Ωsub,1) for near-boundary turbulence and the reddomain (Ωsub,2) where the camera was placed. Particles areseeded to fill the space; within gray areas, particles are seededat higher density to track finer details. . . . . . . . . . . . . . 9Figure 2.5 Overview of our algorithm. At the beginning of each time step,the velocity of each particle is known, the convection-diffusionpart of Navier-Stokes equation is solved by a PSE method(§2.3.1). We then splat the particle velocity to the coarse grid,and make the resulting field divergence-free §2.3.2. The veloc-ity change is interpolated to all particles for divergence correc-tion. Then, for any subdomain Ωsub,i in the field, the correctedparticle momentum is splatted to the grid again and made lo-cally divergence-free. Any particle within a subdomain Ωsub,iis considered a small-scale particle and collects its momentumcorrection from the corresponding domain. . . . . . . . . . . 10Figure 2.6 Impulsively started flow past a sphere at Re=15000 simulatedusing our solver . . . . . . . . . . . . . . . . . . . . . . . . 16Figure 2.7 Vortex shedding in impulsively started flow around a 2m di-ameter sphere at Re=20000. Left: naïve grid-based viscositymodel with h= 0.0625. Middle: naïve solver with h= 0.03125(twice the resolution). Right: our model using a coarse gridof h = 0.125 and a refined grid around the sphere with h =0.03125. The computation time for the left and right simula-tions are roughly equal. . . . . . . . . . . . . . . . . . . . . 17viiiFigure 2.8 For a flow past cylinder simulation at Re = 800, inflow speedU = 5m/s, and cylinder diameter D = 1.6m, the experimentalmodel predicts a Strouhal number of around 0.2 [35], hence ashedding frequency of around 0.625. The above pictures arefrom frame 754 and frame 920 of our simulation with ∆t =0.01, where vortices shed off the same side just reach the endof the domain; this gives a shedding frequency just over 0.6. . 17Figure 2.9 Our method on the same scenario as Fig. 2.7 with differentresolutions. Left: coarse simulation (global h= 0.25, inner h=0.0625). Middle: 2× refinement for both (global h = 0.125,inner h= 0.03125). Right: only the inner grid is refined (globalh = 0.25, inner h = 0.03125). Important fluid features can becheaply captured by increasing just the inner grid resolution. 18Figure 2.10 Rotating fan (each blade is 1.3m long, rotating at one cycleper second) in a constant wind of 2m/s, with ν = 10−4, givingRe ≈ 105. The refined grid has h = 0.005m. Left: zoom-inview at the fan blades. Right: the entire simulation. . . . . . 19Figure 2.11 Smoke rising around a sphere of diameter 0.3m, with ν = 10−4and the inner grid h = 0.005m, Re ≈ 103. Left: simulationwith a uniform FLIP solver. Right: simulation using our solverwith a coarser global grid and a 4× refined grid around thesphere. With approximately the same computation time, oursolver captures structures in the boundary disturbances overthe bottom of the sphere more sharply. . . . . . . . . . . . . 19Figure 3.1 Rising smoke simulations with and without IVOCK (IntegratedVorticity of Convective Kinematics). From top left to bot-tom right,: Stable Fluids, Stable Fluids with IVOCK; BFECC,BFECC with IVOCK; MacCormack, MacCormack with IVOCK;FLIP, FLIP with IVOCK. . . . . . . . . . . . . . . . . . . . 24ixFigure 3.2 Self-advection maps the original velocity field into a rotationalpart and a divergent part, indicated by red and blue arrows re-spectively. Pressure projection removes the blue arrows, leav-ing the rotational part, and illustrating how angular momentumhas already been lost. . . . . . . . . . . . . . . . . . . . . . . 25Figure 3.3 Vorticity confinement (VC) vs. IVOCK. Top row: frame 54of a rising smoke scenario. From left to right, VC parame-ter ε1 = 0.125, ε2 = 0.25, ε3 = 0.5, and SL3-IVOCK. Bottomrow: frame 162. Vorticity confinement tends to create high-frequency noise everywhere, while IVOCK produces a naturaltransition from laminar to turbulent flow with realistic vortexstructures along the way. . . . . . . . . . . . . . . . . . . . . 28Figure 3.4 Vorticity and streamfunction components are stored in a stag-gered fashion on cell edges in 3D (red line), while velocitycomponents are stored on face centers. This permits a naturalcurl finite difference stencil, as indicated by the orange arrows. 31Figure 3.5 Barnes-Hut summation [5] for the boundary values demon-strated in 2D: we construct the monopoles of tree nodes frombottom-up (left), and then for an evaluation voxel (black) in theghost boundary region (grey), we travel down the tree, accu-mulating the far-field influence by the monopoles (blue nodes),and do direct summation only with close enough red cells (right). 32Figure 3.6 Two sequences from 3D rising smoke simulations. Top row:MC-IVOCK with vortex stretching. Bottom row: MC-IVOCKwith vortex stretching switched off. With vortex stretching,vortex rings change their radius under the influence of othervortex rings, these process can easily perturb the shape of vor-tex rings, breaking them to form new vortex structures, whichbrings rich turbulence into the flow field. . . . . . . . . . . . 35Figure 3.7 Left: 2D buoyancy flow simulated with SF. Right: the samewith SF-IVOCK. The resolution and time step were the same;the IVOCK scheme produces a more richly detailed result. . . 36xFigure 3.8 Vortex stretching enhances vortical motion, captured more ac-curately with IVOCK. Transparent renderings illustrate the in-ternal structures. Left and middle-lfet: FLIP. Middle-right andright most: FLIP-IVOCK. . . . . . . . . . . . . . . . . . . . 36Figure 3.9 Left: Frame 150 of an IVOCK simulation with only one multi-grid V-cycle. Right: Frame 150 of an IVOCK simulation wheremultigrid V-cycles are taken reduce the residual by 10−6. . . . 36Figure 3.10 Rising vortex pair initialized by a heat source. Left column: SFwith ∆t = 0.01. Middle column: SF with ∆t = 0.0025. Rightcolumn: SF-IVOCK with ∆t = 0.01. Notice IVOCK preservesvorticity best and produces the highest final position amongthe three. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38Figure 3.11 Vorticity vs. time curve of the 2D vortex pair simulation. WhileIVOCK does not conserve vorticity exactly due to approxima-tions in advection and the streamfunction computation, it stillpreserves significantly more. . . . . . . . . . . . . . . . . . . 38Figure 3.12 IVOCK can be cheaper and higher quality then taking smalltime steps. Left and mid-left: BFECC. Mid-right and right:BFECC-IVOCK with twice as large a time step, computed inmuch less time. . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 3.13 Rising smoke hitting a spherical obstacle: the velocity bound-ary condition is handled entirely by the pressure solve, anddoesn’t enter into the IVOCK scheme at all. This example in-cludes a small additional amount of vorticity confinement toillustrate how the methods can be trivially combined. . . . . . 40Figure 3.14 IVOCK applied to liquids. Top row: dam break simulationsobtained with FLIP (left) and FLIP-IVOCK (right). Bottomrow: moving paddle in a water tank, simulated with FLIP (left)and FLIP-IVOCK (right). In these and several other cases wetested, FLIP-IVOCK is not a significant improvement, presum-ably because interior vorticity is either not present (irrotionalflow) or not visually important. . . . . . . . . . . . . . . . . . 41xiFigure 3.15 Fire simulations making use of our combustion model. Toprow: BFECC. Bottom row: BFECC augmented with IVOCKusing SL3 for vorticity and scalar fields. The temperature isvisualized simply by colouring the smoke according to tem-perature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Figure 4.1 Vortices of raising smoke hitting a moving ball (ball not ren-dered). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Figure 4.2 An overview of our PPPM algorithm, which consists of a far-field construction step and a velocity evaluation step. . . . . . 50Figure 4.3 Relationship between inverse of finite difference operator andGreen’s function. Left: the local inverse matrix defined by A.Middle: the local inverse matrix G constructed using Green’sfunction and the diagonal terms from A. Right: one row of theinverse of G, almost revealed the 7-point stencil. . . . . . . . 55Figure 4.4 Performance of the PPPM fast summation. Computation timegrows linearly with the number of computational elements. . 58Figure 4.5 Accuracy statistics of the PPPM fast summation. . . . . . . . 59Figure 4.6 We switch to a vortex segment representation of vortex parti-cles at the beginning of each time step, move both ends of thevortex segment in the flow, then switch back to vortex blobs. . 62Figure 4.7 Sudden discontinuous motion of vortex segments introducesand amplifies numerical error, which is reduced then by smooth-ing the stretching terms. . . . . . . . . . . . . . . . . . . . . 63Figure 4.8 Vortex shedding process. In our approach, the vorticity strengthis explicitly determined. . . . . . . . . . . . . . . . . . . . . 66Figure 4.9 Rising smoke using different number of vortex particles. Left,2049 particles: middle, 16384 particles; right, 130K particles. 67Figure 4.10 Without the far-field motion, direct summation using a cut-offkernel results in wrong animation. Left: simulation uses cut-off direct summation after 180 time steps. Right: simulationuses PPPM summation after 180 time steps. . . . . . . . . . 68xiiFigure 4.11 Comparison of VIC and PPPM. Top row, sequence of 643 VICsimulation; bottom row: PPPM using the same resolution grid.Notice that the large scale motion of the two simulation matchesbefore the point where turbulent motion becomes dominant. . 68Figure 4.12 Moving objects in slightly viscous flow generate turbulent wakes.Top row: vortex shedding from a moving ball. Bottom row: amoving bunny. . . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 4.13 PPPM Poisson surface reconstruction of: left) a bunny, right)a dragon. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 5.1 3D simulation of flow past sphere at Re=8000. Top: one framefrom the simulated result. Middle: zoom-in of the wake pat-tern. Bottom: near-boundary flow is captured accurately withthe vFLIP solver and seamlessly coupled with the free-spacesolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75xiiiAcknowledgmentsI would like to thank my supervisors Dr. Chen Greif and Dr. Robert Bridson. Theyfed me with interesting questions and provided me helpful insights to my proposedsolutions. I will also thank Dr. Mark Stock: his online summary of vortex methodsliterature has been a great resource for me when I was first exposed to this amazingnumerical method. I also thank Yufeng Zhu, Essex Edwards, and Todd Keelerfor being excellent group-mates, and for the discussions we’ve had to technicalproblems. I also would like to thank all the professors of the IMAGER lab at UBC,for making this lab a great place for computer graphics researches.Thanks also go to my parents for their support; they haven’t been to a collegefor their life and they don’t know science, yet they supported me to pursue myinterests. Thanks go to my friends, too: Yingsai Dong has been my close friend formany years since I started at UBC, we played sports and watched movies together,we celebrated to share happiness and back each other when short for money; heleft Vancouver to work for facebook last year, and I miss him. Thanks go to mygirlfriend, who hasn’t been showing up in my life until recently: you are very verycute.I would also like to thank NSERC (Natural Sciences and Engineering ResearchCouncil of Canada) for funding me to do the published research.xivChapter 1IntroductionBe formless, shapeless, like water. — Bruce Lee.In computer graphics, the animation of fluid such as smoke, fire and wateris usually obtained by numerically advancing fluid quantities with the parametersprovided by artists. To make the numerical techniques a useful tool for users witha minimal mathematical background, unconditionally stable methods [53], despitetheir numerical diffusion, have been widely adopted, trading off the quality of sim-ulations.Fluid motion is simply visually interesting: shear motions in a fluid generateeddies, these eddies merge and split, and move according to each other. These typeof small-scale rotational structures, being highly noticeable in flow (see Fig.1.1),however, are often intractable to capture with the classic semi-Lagrangian fluidsolvers due to numerical diffusion.To combat the numerical diffusion problem, many methods have been proposedeven just within computer graphics, such as higher order advection solvers [31,50], momentum-conserving advection [34], and the non-diffusive Fluid ImplicitParticle (FLIP) method [68].As an alternative to these velocity-pressure schemes, researchers in and out-side graphics field have been using vortex methods to track the turbulent motion offluids. Vorticity is highly related to the swirling motion of flows; tracking the evo-lution of vorticity directly can often better generate and preserve detailed fluid fea-1Figure 1.1: The chaotic motion of eddies produced by the storms in the at-mosphere of Jupiter. Courtesy NASA/JPL-Caltech.tures, with less computational elements than a typical computer graphics velocity-pressure fluid solver might require.An early milestone in computer graphics used the vortex method and the Vortex-in-Cell (VIC) computation to reproduce the turbulent eddies on Jupiter [65]. Later,Park and Kim demonstrated that complex fluid motions can be simulated withonly 8K vortex particles [44]. Vortex methods have also been used to study highReynolds number flows [12, 55, 57].However, many disadvantages of the vortex methods have also been observed.Calculating the flow velocity from vorticity requires O(N2) computations for Nvortex elements with a straightforward application of the Biot-Savart law. Thenumber of vortex elements might grow exponentially during an entire simulation,making the computation intractable even when O(N) fast summation is used. Lastbut not least, boundary conditions in vortex methods require non-trivial treatment,and as a result, very limited boundary conditions have been applied in classicalvortex method simulations.We observe that the creation, evolution and interaction of vorticity elements is2highly related to the turbulent, swirling feature of flow motions, but certain trade-offs have to be made to the vortex method for its efficient, flexible and robust appli-cation in computer graphics. This thesis discusses the success we have achieved inadapting the vortex method regarding the three aspects, and points the way forwardto further improvements.3Chapter 2Resolving Fluid Boundary Layerswith Particle Strength Exchangeand Weak Adaptivity2.1 IntroductionIn computer graphics, the incompressible Navier-Stokes equations are often usedto produce realistic fluid animation. Storing fluid quantities on a Cartesian gridand performing pressure projection to handle incompressibility and boundary con-ditions, Eulerian approaches have proved their merit in scalably handling complexboundary shapes and maintaining stability with large time steps [7].Despite the ease of use and implementation of Eulerian approaches, their abil-ity to simulate high Reynolds flow remains a problem. In this type of flow, fluidmotion is often strongly influenced by boundary layer dynamics. A boundarylayer of vanishing thickness usually cannot be resolved at practical grid resolu-tions, making the no-stick or no-slip boundary conditions both diverge from phys-ical predictability. Besides the possibility of producing inconsistent results undergrid refinement, the flow motion is poorly or not at all controlled by changing theReynolds number, as shown in Fig. 2.2: simulations at low resolution, which don’tresolve the boundary layer, cannot hope to resemble the results produced by high4Figure 2.1: Simulation of a rotating fan in wind, Re ≈ 105, visualized withsmoke marker particles. Large scale motion and detailed boundary layerdynamics are cheaply and easily coupled using our method, producinghigh quality fluid animations at low cost. Top: top view of fan wake atframes 80, 120 and 200. Bottom: side view at frame 300.resolution simulations.Fluid-structure interactions in slightly viscous flow has been successfully mod-eled by vortex methods [9, 56]. Boundaries are viewed as generators of vorticity invortex methods. However, solving the boundary integral equation for general ge-ometry in 3D is non-trivial, and concerns also remain about how to reliably achievestability with 3D vortex stretching [20, 66]: we have found it hard in practice togenerally adopt vortex methods for computer graphics.On the other hand, while the variational framework [6] with FLIP advection[68] is capable of simulating free-slip boundaries nicely, with insufficient grid res-olution the momentum exchange near (more physical) no-slip boundaries is onlypoorly resolved.To capture near-boundary flow accurately without incurring the same expensefor the domain as a whole, adaptive methods have been proposed for graphics bymany researchers, e.g. using adaptive grids [2, 37, 51] or domain decomposition[21]. While the former category has advantages in more smoothly varying thedegree of refinement, there is extra overhead in mesh generation and memory in-dexing. On the other hand, domain decomposition solvers require extra iterationsper time-step to couple the solutions between different domains.5Figure 2.2: Vortex shedding from a cylinder in 2D at Re=15000, zoom-inview near the boundary. Left column: simulation results at low resolu-tion (approximately 50× 50 grid cells shown). Middle column: simu-lation obtained with 4× resolution. Right column: simulation obtainedwith our method, with the same coarse grid but a 4× refined grid over-laid just around the solid. Our method produces results visually consis-tent with the high resolution reference because only the boundary layerneeds that resolving power.We take inspiration from how, in FLIP, the particles carry flow velocities andmove with the flow, while an Eulerian grid is used to adjust the particle velocitiesto a divergence-free state. In this paper, we extend this philosophy to arrive at anefficient, easy to implement, highly adaptive fluid solver:• The classic FLIP scheme is augmented with a particle strength exchange(PSE) method [39] to solve the convection-diffusion part of the Navier-Stokes equations.• By seeding extra ghost boundary particles at the boundary, our solver cap-tures the boundary layer dynamics more accurately, producing results visu-ally consistent with higher resolution simulations (cf. Fig. 2.2).• Our time integration weakly couples the regionally refined solutions effi-ciently, alleviating the need for sophisticated global solvers.• The regional refinements can be placed arbitrarily, even overlapping each6Table 2.1: Symbol abbreviations used throughout this paper.ν Kinematic viscosity of the fluidup,i Velocity stored on particle ixi Position of particle iug(x) Velocity sampled on grid at position xΩ Fluid domainΩsub,i ith Sub-domain with particle-mesh refinementWh(x) Kernel used to spread a particle quantity to the gridHτ(x) Heat kernel for particle momentum exchangeh Local grid cell spacingFigure 2.3: Our method captures the velocity field near the boundary effi-ciently and with high apparent fidelity. Left: zoom-in near the boundaryflow. Right: the entire simulation, Re=15000.other, without any geometric operations to merge them together, simplifyingmesh generation for spatially adaptive solvers.2.2 Related workBoundary layers are often the source of turbulence in high Reynolds number flows.Pfaff et al. [45] attempted to capture this effect by seeding vortex particles from aprecomputed artificial boundary layer to enhance a coarse simulation. However,their approach is limited to static boundaries and may be less physically plausiblethan desired as the vortex particles are seeded by chance. In contrast, our methodaims to solve the viscous dynamics at the right scale, alleviating the limitation tostatic boundaries and improving direct physical control of the result.When modelling fluid motion with vortex elements, vorticity can be seedednear the boundary by diffusing the boundary vortex sheet [44, 56]. In these ap-proaches, the unknown vortex sheet strength is determined by solving boundary7integral equations. Besides the expense of the GMRES solve for the vortex streetstrength, there are questions of solvability of the equations themselves for arbi-trary topology under general motion, which is significant for computer graphicsapplications.Particles, along with Eulerian velocity projections, have been widely used incomputer graphics to capture fluid features which may fall between grid samplesand be smoothed away by purely Eulerian solvers, e.g. FLIP [68], or derivative par-ticles [52] and the closely related APIC [29]. We further extend this concept: in ourmethod, particles are used to represent the change of flow momentum due to diffu-sion and external forces, while multiple-resolution Eulerian grids bounding differ-ent subregions are used to project particle momentum towards their divergence-freestate (along with a global coarse grid).Our scheme is related to the Iterated Orthogonal Projection (IOP) frameworkproposed by Molemaker et al. [41] but differs in two aspects: instead of re-projectingthe field solution globally in consecutive iterations, our projection is only appliedto a set of locally refined sub-domains; during each projection our pressure dis-cretization respects the solid boundary conditions as per Batty et al. [6], while inIOP, boundary conditions are only imposed in a separate step of the iteration.Chimera grids proposed by English et al. [15] are a promising tool for adap-tive, large-scale fluid computations, but require nontrivial domain discretization tomerge grids together, and a relatively expensive global solver. While the globalsolve may be critical for water simulations, in this paper we focus our attention onpurely gaseous phenomena and demonstrate a faster and simpler weak coupling ofgrids is effective.Domain decomposition schemes have been used as preconditioners for iterativelinear solvers within pressure projection (e.g. [13]), or as a useful tool to decouplefluid features with different solution methods [21]. Our method also makes useof domain decomposition concepts: we solve the convection-diffusion and forcecalculation part on particles, and use multiple Eulerian grids as an auxiliary toolto perform velocity projection to obtain an adaptive refined and weakly coupledsolution in the fluid domain.Stock et al. [57] proposed an efficient one-way coupling to simulate rotor wake.The fluid state is first updated with a vortex particle time step, which provides the8ΩΩ ,2Ω ,1Figure 2.4: A typical domain construction used in our method. Ω indicatesthe global Eulerian domain where fluid motion is loosely captured. Eu-lerian subdomains with finer resolution can be placed anywhere to en-hance the simulation quality locally, such as the green domain(Ωsub,1)for near-boundary turbulence and the red domain (Ωsub,2) where thecamera was placed. Particles are seeded to fill the space; within grayareas, particles are seeded at higher density to track finer details.velocity boundary condition for the Eulerian sub-domain to advance itself. Vor-tex particles in the Eulerian sub-domain can then interpolate vorticity changes forthe next time-step. In contrast, our space-filling Lagrangian particles carry fluidmomentum instead of vorticity, to avoid the potential complexities associated withvortex stretching and viscous boundary conditions in general three-dimensionalflows.Multigrid solvers are becoming increasingly popular in graphics to acceleratethe pressure solve (e.g. [18, 40], [3]) We likewise use multigrid for the variousgrids in our solver.While Lentine et al. [33] used multi-level grids to obtain detailed fluid anima-tions efficiently, their algorithm is limited to purely Eulerian schemes, whereas oursolver is designed for hybrid particle-mesh methods like FLIP, which offer lessnumerical diffusion.A common technique in industry is to inject additional detail with vorticityconfinement [16, 54]) or post-process turbulence synthesis (e.g. [32, 48]). We didnot use such methods, focusing on a more direct approach of better solving the9Figure 2.5: Overview of our algorithm. At the beginning of each time step,the velocity of each particle is known, the convection-diffusion part ofNavier-Stokes equation is solved by a PSE method (§2.3.1). We thensplat the particle velocity to the coarse grid, and make the resulting fielddivergence-free §2.3.2. The velocity change is interpolated to all par-ticles for divergence correction. Then, for any subdomain Ωsub,i in thefield, the corrected particle momentum is splatted to the grid again andmade locally divergence-free. Any particle within a subdomain Ωsub,i isconsidered a small-scale particle and collects its momentum correctionfrom the corresponding domain.Navier-Stokes equations at higher Reynolds numbers, where the look can be con-trolled primarily with actual physical parameters. However, procedural techniqueslike these could be added as an additional layer to our method to produce evengreater details.2.3 VFLIP and weakly higher resolution regionalprojection (WHIRP).The incompressible Navier-Stokes equations we approximately solve are:∂u∂ t+(u ·∇)u =− 1ρ∇p+ν∆u+ f ,∇ ·u = 0.(2.1)We adopt the usual FLIP framework to handle the material derivative Du/Dt onthe left hand side, storing velocity on particles which are moved in a grid-basedvelocity field. The pressure gradient and incompressibility condition (togetherwith boundary conditions) are handled by a separate pressure projection step, aug-mented in this paper with multiple grids (see figure 2.4). We use particle strengthexchange (PSE) for the viscous term, and integrate body forces f with an Euler10step split from the rest of the time integration as usual.An overview of our algorithm is given in Fig. 2.5, while each subroutine calledin a time step is listed in Alg. 1. Technical details of each subroutine are given inthe corresponding subsections, §2.3.1, §2.3.2,§2.3.3 and §2.3.4.Algorithm 1 TimeStep(∆t )1: // Convection-diffusion §2.3.12: For each particle i3: up,i = PSE(∆tν);4: xi = ForwardTrace(ug, ∆t, xi)5: // Hierarchical projection §2.3.26: In fluid domain Ω7: particles = Collect particles in Ω8: DivergenceFreeUpdate(particles)9: For each sub-domain Ωsub,i in dx descending order10: particles = Collect particles in Ωsub,i11: DivergenceFreeUpdate(particles)12: SeedAndDeleteParticles // §2.3.4Algorithm 2 DivergenceFreeUpdate(particles)1: u∗g = Splat particle velocities to grid cells2: ug = Project(u∗g)3: δu = ug−u∗g4: For each particle i5: up,i = up,i+δu(xi)2.3.1 Solving the convection-diffusion equation with ghost particlesGiven particles and their velocity at time step n, the convection-diffusion part ofthe Navier-Stokes equation,DuDt= ν∆u, (2.2)can be solved efficiently with a so-called Particle Strength Exchange method. Todeal with boundary objects, at each time step we seed ghost particles randomly in a2h band inside the boundary (where h is the spacing of the grid covering the solid,and with the same density as regular FLIP particle seeding), and set their velocity11to that of the solid object. Given particle velocity up,i, the updated particle strengthis then calculated using the stable formulaun+1p,i = unp,i+∑ j∈η(up, j−up,i)Hτ(xi− x j)∑ j∈η Hτ(xi− x j)(2.3)where η is the neighborhood in which we look for particles (a box of size 2h), andHτ(x) is the heat kernel with τ = νδ t, which readsHνδ t(x) =1(4piνδ t)d/2exp(−‖x‖24νδ t)(2.4)The (4piνδ t)−d/2 factor can be elided due to the normalization in Equation 2.3.Once the particle velocites have been diffused, particles are passively advectedthrough the divergence-free flow field. We used Ralston’s 3rd order Runge-Kuttatime integrator [47], as the most efficient optimal RK scheme with a stability regionoverlapping the pure imaginary axis (for rotations):u′i = SampleVelocity(xi)u′′i = SampleVelocity(xi+0.5∆tu′i)u′′′i = SampleVelocity(xi+0.75∆tu′′i )xn+1i = xni +29∆tu′+39∆tu′′+49∆tu′′′(2.5)When sampling velocity from the velocity buffer, we use trilinear interpolation.In the case where multiple overlapping grids are detected, we take the grid withthe finest resolution. While this sounds overly simple, and potentially introducesdiscontinuities at the edges of grids, we have not observed any noticeable artifactsas typically the velocity jumps would be small, and it is only the integral in time ofthe velocity (the particle trajectories) which we observe. It would not be difficult touse a partition-of-unity blend to get a smoother velocity at grid edges if for somereason this was necessary. After advection, we further update the particle velocitieswith forcing terms such as buoyancy before pressure projection.122.3.2 Regional projection for particle velocity correctionOnce the post advection particle velocities are known, we use a regional projectionmethod to correct the velocity field to an approximately divergence-free state. Fol-lowing Batty et al. [6], finding the inviscid divergence-free projection of an inputvelocity field u∗ is equivalent to a minimization problem:minpˆΩ12ρ‖u∗− ∆tρ∇p‖2−ˆΩ12ρ‖u∗‖2+ˆSpnˆ ·∆tupre, (2.6)where upre is typically a prescribed velocity at solid or fluid domain boundaries S.As outlined in Alg. 2, in our regional projection, we first splat all particle veloc-ities to the Eulerian domain Ω with a wide kernel Wh(x), as if we were computingthe large eddy filtered version of the velocity field:ug(x) =∑ j u jWh(x− x j)∑ j Wh(x− x j). (2.7)Here x is a fixed grid cell position, and x j and u j are the the position and velocityof particle j respectively.In our implementation, the usual bi-/trilinear hat function kernel was used:Wh(x) =d∏θ=1φ(xθh),withφ (r) =1−|r|, r ∈ [−1,1) ,0, else,(2.8)where d = 2,3 is the spatial dimension, h is the kernel width which simply equalsthe grid spacing, and xθ is the θ th component of vector x.After we splat the particle velocity to the grid, the velocity field on the gridis made divergence-free by pressure projection. Our projection uses a standardstaggered grid and face-weighted pressure discretization similar to Batty et al. [6].The velocity change is then interpolated to the particles to correct their velocities.At this stage the velocity defined on the particles is (approximately) divergence-free at the coarse scale. To increase the flow fidelity in regions of interest (nearboundaries or cameras, for example), further projections with higher resolution13in each subdomain are performed with essentially the same procedure. We splatthe coarse-corrected particle velocity to the corresponding refined domain usingthe kernel Wh(x) with the smaller refined h, make the grid-based velocity fielddivergence-free, and interpolate the velocity change back to the particles in thatregion. In each subdomain projection, at the domain boundaries we restrict thenormal component of velocity not to change (as established by splatting from theparticles). We repeat this process until all the regional refined domains are pro-cessed.During each projection, the Poisson equation for pressure is solved with amulti-grid preconditioned conjugate gradient method, whose implementation de-tails are given in §2.3.3.2.3.3 Our pressure solverSimilar to McAdams et al. [40] and Setaluri et al. [51], we use a multigrid precon-ditioned Krylov solver for our pressure Poisson equation. Our multigrid solver isconstructed semi-algebraically: we construct RL (the restriction matrix mappingthe residual of level L to a coarser level) and PL (the prolongation matrix usedto add coarse level corrections back to the level L solution) using geometric infor-mation given by the fluid domain, but use matrix multiplications to construct thecoefficient matrix AL+1 of each coarsening level L+1 (Alg. 3).Algorithm 3 MultiGridCoarsening(AL,PL,RL, AL+1)1: for each cell in coarse level2: i← index of this cell3: for each sub cell in fine level4: j← index of this sub cell5: if AL( j, j) 6= 06: RL(i, j) = 1/2d //d = 2,3 dimension of the problem7: PL( j, i) = 1.08: AL+1← 0.5RLALPLAs shown in Alg. 3, our restriction and prolongation matrices are transposesof each other, up to a scale factor, with restriction equivalent to simply averagingresiduals across the fine level grid cells covered by a coarse grid cell, and prolon-gation equivalent to taking the nearest neighbour coarse value (i.e. ”parent” value);14this is a ”aggregation”-style multigrid with constant weights. For pure multigrid,this does not produce optimal convergence; we improve the scheme by applying ascaling of 12 to the coarse level coefficient matrix, which, away from solid bound-aries, gives exactly the same coefficient matrix as if we were to rediscretize thePDE on the coarse grid. For cases with solid boundaries, this coarsening strat-egy acts as a volume weighted rediscretization at the coarse level. Our coarseningstrategy doesn’t increase the size of the discrete Poisson stencil: a coarse cell has atmost six neighbors with non-zero coefficients. This in turn allows for the usual red-black ordering to parallelize Gauss-Seidel sweeps, and in general makes smoothingfar more efficient. For other implementation details of multigrid-preconditionedKrylov solvers, please refer to McAdams et al. [40].In all our experiments, the multigrid-preconditioned conjugate gradient solverappeared to converge with a constant rate of approximately 0.1 per iteration, asshown in Table 2.2, though more complex schemes may be necessary for extremelycomplex solid geometry.Table 2.2: The Semi-AMGPCG solver applied to a Poisson problem with sin-gle impulse, in a unit cube domain with a zero Neumann boundary con-dition along a sphere with diameter 0.2 in the center, and zero Dirichletboundary conditions at the edges of the cube. The number of iterationsrequired to reduce the residual by a factor of 10−10 is essentially indepen-dent of the grid size, providing one extra digit of accuracy per iteration.Convergence rate of the AMGPCG solverResolution 643 963 1283 1923 2563Convergence criterion 10−10||r.h.s||∞# of iterations 10 11 11 11 112.3.4 Seeding and deleting particlesAt the end of each time step, particles are seeded and deleted as necessary to main-tain reasonable sampling of the fluid in each grid cell, according to the highestresolution grid nearby. Newly seeded particles take their velocity by interpolationfrom the finest available grid. We delete particles randomly from cells where thecount exceeds a threshold (16 in our examples), and seed new particles where the15Figure 2.6: Impulsively started flow past a sphere at Re=15000 simulated us-ing our solvercount is below another threshold (8 in our examples).Referring to figure 2.4, the grey areas surrounding (but not inside) refined sub-domains are treated with the finer grid spacing: we maintain a higher particle den-sity in the grey areas as there is a chance that the fluid in the grey region will beadvected into the refined subdomain during the course of the next time step. Forour examples we used a bandwidth of 3h, where h is the refined grid spacing, forthis intermediate zone but it could of course be determined more dynamically andfrugally based on the local velocity and the time step size.As an alternative to avoid higher particle counts, we have experimented withonly seeding extra particles in the refined regions, not the grey neighborhoods, butwidening the support of the particle-to-grid kernel near the edge of the refined gridto avoid gaps in the data. This exchanges more work for less storage, which maybe preferable in some cases, though the results are similar.2.4 Results and discussionWe parallelized the proposed methods on a desktop machine with an Intel(R)Core(TM) i7-3930K CPU and 32 GB RAM. We compared our new solver to a typ-ical inviscid FLIP solver as used in graphics, as well as a “naïve” viscous Navier-Stokes solver where viscosity is added with an explicit-in-time finite difference(since we are interested in very high Reynolds numbers, this explicit step was al-ways stable). Simulations were visualized by judiciously injecting smoke markerparticles which are passively advected with the grid velocity field. The accompa-nying video shows comparisons between the methods for different resolutions anddifferent Reynolds numbers in a variety of examples.In Fig. 2.6, we show a flow simulation frame obtained with our solver at16Figure 2.7: Vortex shedding in impulsively started flow around a 2m diametersphere at Re=20000. Left: naïve grid-based viscosity model with h =0.0625. Middle: naïve solver with h = 0.03125 (twice the resolution).Right: our model using a coarse grid of h = 0.125 and a refined gridaround the sphere with h = 0.03125. The computation time for the leftand right simulations are roughly equal.Figure 2.8: For a flow past cylinder simulation at Re= 800, inflow speed U=5m/s, and cylinder diameter D= 1.6m, the experimental model predictsa Strouhal number of around 0.2 [35], hence a shedding frequency ofaround 0.625. The above pictures are from frame 754 and frame 920of our simulation with ∆t = 0.01, where vortices shed off the same sidejust reach the end of the domain; this gives a shedding frequency justover 0.6.Reynolds number of 15000: the flow separates from the boundary at small an-gle, it remains laminar with structured vortices for about one diameter, and thenbecomes turbulent. This is similar to figure 55 on page 34 of An Album of FluidMotion [61].As shown in Fig. 2.7, with the naïve uniform-resolution grid-based viscositymodel, a much higher resolution is necessary to capture important fluid motionssuch as vorticity detaching from the boundary layer. With our sub-grid particle-based viscosity model and regional refinement, the characteristic vortex sheddingcan be achieved at much lower cost.In Fig. 2.8 we performed a 2D flow past cylinder simulation. The cylinderis only about 12 cells wide in the coarse grid; a 4× local refinement is appliedto resolve the boundary layer. The shedding frequency produced by our solver17Table 2.3: Vortex shedding frequency statistics from numerical simulationsVortex shedding frequency fν ofa flow-past-sphere simulation in 2D with various Strouhal numbersReynolds number 400 1600 16,000Strouhal number 0.188 0.195 0.197Ground truth fν 0.117 0.122 0.62Numerical fν 0.114 0.121 0.611Figure 2.9: Our method on the same scenario as Fig. 2.7 with different res-olutions. Left: coarse simulation (global h = 0.25, inner h = 0.0625).Middle: 2× refinement for both (global h = 0.125, inner h = 0.03125).Right: only the inner grid is refined (global h = 0.25, inner h =0.03125). Important fluid features can be cheaply captured by increas-ing just the inner grid resolution.matches the experimental model. We also studied this type of flow with variousStrouhal numbers: results are given in Table 2.3.In Fig. 2.9 we take the same scenario (impulsively started flow around a sphereat Re = 20,000) but look at changing global and inner grid resolution in ourmethod. Vortex shedding is largely dominated by how well a solid boundary isresolved. Our method can resolve some sub-grid scale viscosity effects, but it failsto predict the near-boundary behavior when the refined grid is too coarse to resolveit. However, increasing the resolution of the local grid just around the boundaryimproves the simulation quality greatly without requiring refinement in the rest ofthe domain.When the resolution of the coarse domain isn’t enough to represent the solidobject, our regional refinement may still be enough to capture the solid boundarywell (Fig. 2.1, Fig. 2.10). Vortex shedding from solid objects is essential to fluidanimation in these scenarios, which can be faithfully captured by our method.Our method makes it much cheaper to adequately resolve the viscous bound-ary layer, locally and on-the-fly, alleviating the limitation to parametrize boundary18Figure 2.10: Rotating fan (each blade is 1.3m long, rotating at one cycle persecond) in a constant wind of 2m/s, with ν = 10−4, giving Re ≈ 105.The refined grid has h= 0.005m. Left: zoom-in view at the fan blades.Right: the entire simulation.Figure 2.11: Smoke rising around a sphere of diameter 0.3m, with ν = 10−4and the inner grid h = 0.005m, Re≈ 103. Left: simulation with a uni-form FLIP solver. Right: simulation using our solver with a coarserglobal grid and a 4× refined grid around the sphere. With approxi-mately the same computation time, our solver captures structures inthe boundary disturbances over the bottom of the sphere more sharply.layer with far field flow motions as in Pfaff et al. [45]. It is thus suitable for moregeneral simulation scenarios, e.g. buoyancy-driven plumes as shown in Fig. 2.11.In Table 2.5, we notice that for the inner grid of 2003, the simulation time isvery close while the outer grid size varies from 100x50x50 to 150x75x75. This isdue to the reason that in our experiments, as the resolution of inner grid increases,the particles assigned for the inner domain increases, as a result, the computationtime became dominated by the viscous diffusion part, where each particle’s veloc-ity is adjusted by the kernel summation. In practice, diffusion only needs to beresolved up to a narrow band near the boundary, hence, the cost function can bequantified asT = c1Ng+ c2Nb (2.9)where, Ng is the total number of grid cells in the entire domain, including the coarse19Table 2.4: Timings of various simulations associated with this paper, mea-sured on a desktop machine with Intel(R) Core(TM) i7-3930K CPU and32 GB RAM.Reynolds Simulation TimeExample Method ∆t Outer Grid Size Inner Grid Size Number (seconds per frame)Flow past sphereNaïve solver0.03 0.0625 N/A 1,000 12.520.03 0.03125 N/A 20,000 73.46Our method0.03 0.125 0.0625 1,000 11.210.03 0.0625 0.03125 20,000 70.760.03 0.125 0.03125 15,000 56.07Rotating fan in windOur method0.03 0.041 0.02 ≈ 105 120.40Rising smoke plume 0.05 0.1 0.03 ≈ 103 79.10Flow past sphereFLIP 0.03 0.0625 N/A N/A 11.31FLIP + WHIRP 0.03 0.125 0.0625 N/A 4.06VFLIP 0.03 0.0625 N/A N/A 15.41VFLIP + WHIRP 0.03 0.125 0.0625 N/A 4.20Table 2.5: Timings of flow past sphere simulations with only differences inouter and inner grid sizes, measured on a desktop machine with Intel(R)Core(TM) i7-3930K CPU and 32 GB RAM.Simulation Time(seconds per frame)Outer Grid Size100x50x50 150x75x75 200x100x100InnerGridSize66x66x66 12.19 19.71 34.36100x100x100 28.99 35.46 50.52200x200x200 166.14 168.64 185.57and subdomains, as the fluid advection and projection operations shall take similarcost for every grid cell. On the other hand, the diffusion scales linearly with thenumber of voxels in a narrow band near the boundary. We fitted our experimentaldata to get c1 = 0.0000126 and c2 = 0.000120. This shows that the particle diffu-sion time becomes dominant for highly refined simulations. We believe that ratherthan the current CPU parallelizeed code, a GPU parallelized code could be usedto significantly reduce the computation time for the particle diffusion, as we did inour 2D examples.2.5 Conclusions, limitations, and future workWe extended incompressible FLIP in two ways, with a particle-based viscositymodel that can resolve momentum exchange at higher resolution than the grid, andwith an extremely simple regional projection method (together with particle seed-20ing/deletion rules) to adaptively refine in important areas, without needing morecomplex grid or mesh structures to globally couple different resolutions. Very littleextra code beyond a standard FLIP simulator is required to implement this paper.Moreover we showed that matching the qualitative look of some high Reynoldsnumbers scenarios hinges just on resolving the viscous boundary layer, even whilethe rest of the domain can use a coarse grid.The typical inviscid solver used in graphics doesn’t take viscosity into accountat all, and for example can only lead to vortex shedding through numerical errorsrelated to grid size and time step. Both the naïve method and our method will failto give characteristic results when the boundary layer isn’t resolved, again leadingto simulations whose look is more controlled by numerical errors related to gridsize and time step than physical parameters. However, our method allows for re-solving boundary layers much more efficiently, at which point physical parametersare at last a useful control and further refinement will reliably give similar but moredetailed results, as artists generally hope for.Our solver is still limited formally to first order accuracy in time, and the PSEapproach is likely only first order in space: while we feel it’s valuable for coarsegrids and large time steps, it is not an efficient approach for convergence to a fullyaccurate solution. The PSE part is also inappropriate for highly viscous flows,where an implicit grid-based solver is almost certainly the better option.The global time step is also a limitation of this approach to spatial-only adaptiv-ity. The time step selected in graphics simulations, when multiple steps per frameare taken, is frequently determined by restricting the CFL number, i.e. limiting thenumber of grid cells a fluid particle may travel through in one time step. In refinedsubdomains with smaller grid cells, this is a much more stringent restriction thanis needed for the coarse global grid, wasting some computational effort.We have not considered fluid phenomena beyond simple smoke scenarios. Fireand explosions, where nonzero divergence may be prescribed by the pressure pro-jection, would require some adjustment to our regional projection scheme. Freesurface liquids, such as commonly used for water, would require even more thoughtif regional grids overlap the free surface.There are several other avenues for further research. The convection-diffusionpart we currently handle with PSE might be improved with other ideas from grid-21free methods; we could also look at using an effective eddy viscosity to bettermodel unresolved turbulence. Our regional projection is not limited to regulargrids, but could also use unstructured tetrahedral meshes to conform more accu-rately to solid boundaries (while avoiding the need for a global, conforming, tetra-hedral mesh, which could be a useful savings for simulations involving movingrigid bodies). Last but not least, it is possible to also track the vorticity on eachparticle to correct the missing angular momentum during particle advection, simi-lar to Zhang et al. [67]; replacing the FLIP framework with APIC [29] might alsoallow for better tracking of vorticity.22Chapter 3Restoring the Missing Vorticity inAdvection Projection FluidSolversTable 3.1: Algorithm abbreviations used through out this paper.IVOCK The computational routine (Alg.4) correctingvorticity for advectionSF Classic Stable Fluids advection [53]SF-IVOCK IVOCK with SF advectionSL3 Semi-Lagrangian with RK3 path tracingand clamped cubic interpolationBFECC Kim et al.’s scheme [31], withextrema clamping([50])BFECC-IVOCK IVOCK with BFECC advectionMC Selle et al.’s MacCormack method [50]MC-IVOCK IVOCK with MacCormackFLIP Zhu and Bridson’s incompressible variantof FLIP [68]FLIP-IVOCK FLIP advection of velocity and density,SL3 for vorticity in IVOCK.23Figure 3.1: Rising smoke simulations with and without IVOCK (IntegratedVorticity of Convective Kinematics). From top left to bottom right,:Stable Fluids, Stable Fluids with IVOCK; BFECC, BFECC withIVOCK; MacCormack, MacCormack with IVOCK; FLIP, FLIP withIVOCK.3.1 IntroductionIn computer graphics, incompressible fluid dynamics are often solved with fluidvariables stored on a fixed Cartesian grid, known as an Eulerian approach [53].The advantages of pressure projection on a regular grid and the ease of treatingsmooth quantities undergoing large deformations help explain its success in a widerange of phenomena, such as liquids [19], smoke [16] and fire [43]. Bridson’s textprovides a good background [7].In Eulerian simulations, the fluid state is often solved with a time splittingmethod: given the incompressible Euler equation,∂u∂ t+(u ·∇)u =− 1ρ∇p+ f∇ ·u = 0(3.1)fluid states are advanced by self-advection and incompressible projection. First24Original velocity fieldri i l l it fi l Velocity field after advectionl it fi l ft r tiVelocity field after pressure projectionl it fi l ft r r r r j tiFigure 3.2: Self-advection maps the original velocity field into a rotationalpart and a divergent part, indicated by red and blue arrows respectively.Pressure projection removes the blue arrows, leaving the rotational part,and illustrating how angular momentum has already been lost.one solves the advection equationDuDt= 0 (3.2)to obtain an intermediate velocity field u˜, which is then projected to be divergence-free by subtracting ∇p, derived by a Poisson solve from u˜ itself. We formallydenote this as un+1 = Proj(u˜).Self-advection of velocity, without continuous application of the pressure gra-dient as is typical of the first-order time-splitting schemes used in graphics, dis-regards the divergence-free constraint, allowing some of the kinetic energy of theflow to be transferred into divergent modes which are lost in pressure projection.We focus in particular on rotational motion: self-advection can sometimes cause anoticeable violation of conservation of angular momentum, as illustrated in Fig.3.2.Despite many published solutions in improving the accuracy of advection schemein an Eulerian framework (e.g. [34, 50]), or reducing the numerical diffusion withhybrid particle-in-cell solvers [68], even with higher-order integration [13], only afew papers have addressed this particular type of numerical dissipation.Retaining the velocity-pressure approach but changing the time integration toa symplectic scheme, Mullen et al. [42] were able to discretely conserve energywhen simulating fluids. However, the expense of their non-linear solver for Crank-Nicolson-style integration, and the potential for spurious oscillations arising fromnon-upwinded advection, raise questions about its practicality in general.Another branch of fluid solvers take the vorticity-velocity form of the Navier-25Stokes equation, exploiting the fact that the velocity field induced from vorticity isalways divergence free. These solvers advect the vorticity instead of the velocity ofthe fluid, preserving circulation during the simulation. Vorticity is usually trackedby Lagrangian elements such as vortex particles [44], vortex filaments [64] or vor-tex sheets [8]. Apart from the fact that the computational cost of these methodscan be a lot more expensive per-time-step than a grid-based solver (three accu-rate Poisson solves are required instead of just one), there remain difficulties inhandling solid boundaries and free surfaces (for liquids), and users may find itless intuitive working with vorticity rather than velocity in crafting controls for artdirection. However, researchers have observed that vortex methods can better cap-ture important visual details when running with the same time step and advectionscheme [65].In this paper, we combine our observation from vortex methods with Euleriangrid-based velocity-pressure solvers to arrive at a novel vorticity error compensa-tion scheme. Our contributions include:• A new scheme we dub IVOCK (Integrated Vorticity of Convective Kinemat-ics), that tracks and approximately restores the dissipated vorticity duringadvection, independent of the advection method or other fluid dynamics inplay (boundary conditions, forcing terms) being used.• A set of novel techniques to store and advance vortex dynamics on a fixedspatial grid, maximizing the accuracy while minimizing the computationaleffort of IVOCK.• Upgraded classic fluid solvers with IVOCK scheme, documenting the com-bination of IVOCK with different advection methods and different types offluids.• A new vortex stretching model in a non-divergence-free environment forcompressible vortex flows occuring in volumetric combustion.3.2 Related workStam’s seminal Stable Fluids method [53] introduced a backward velocity trac-ing scheme to solve advection equation, making physically based fluid animation26highly practical for computer graphics. With its unconditional stability, tradingaccuracy for larger time step, it provided a basis for most subsequent grid-basedfluid solvers in graphics, for fluid phenomena such as smoke [16], water [19], thinflames [43], and volumetric combustion [17].However, the first order accuracy in both time and space of Stable Fluids man-ifests in strong numerical diffusion/dissipation. Many researchers have taken it asa basic routine to build higher-order advection schemes. Kim et al. [31] proposedthe BFECC scheme to achieve higher-order approximation by advecting the fieldback and forth, measuring and correcting errors. Selle et al. [50] eliminated thelast advection step of BFECC to arrive at a cheaper MacCormack-type method,and also introduced extrema clamping in BFECC and MacCormack to attain un-conditional stability, at the cost of discontinuities in velocity which can sometimescause visual artifacts. In this paper, we use these advection schemes as examplesto which we apply IVOCK, and show improvements with the IVOCK scheme forsmoke and volumetric combustion simulations.Hybrid particle-grid methods have been introduced to further reduce numer-ical diffusion in adection, notably Zhu and Bridson’s adaptation of FLIP to in-compressible flow [68]. Although FLIP is almost non-diffusive for advection, thevelocity-pressure solver nevertheless may dissipate rotational motion as shown infigure 3.2, regardless of the accuracy of advection. FLIP cannot address this issue;in this paper we show our FLIP-IVOCK scheme outperforms FLIP in enhancingrotational motions for smoke.Resolution plays an important role in fluid animations. McAdams et al. [40]proposed multigrid methods to efficiently solve the Poisson equation at a uniformhigh resolution. Losasso et al. [37] introduced a fluid solver running on an octreedata structure to adaptively increase resolution where desired, while Setaluri etal. [51] adopted a sparse data approach. In our work, multigrid is employed for thevorticity-stream function solver while an octree code is used to impose domain-boundary values using Barnes-Hut summation.Vortex dynamics has proved a powerful approach to simulating turbulence. La-grangian vortex elements such as vortex particles [44], vortex filaments [64], orvortex sheets ([8], [46]) have been used to effectively model the underlying vor-ticity field. Recently, Zhang and Bridson [66] proposed a fast summation method27Figure 3.3: Vorticity confinement (VC) vs. IVOCK. Top row: frame 54 ofa rising smoke scenario. From left to right, VC parameter ε1 = 0.125,ε2 = 0.25, ε3 = 0.5, and SL3-IVOCK. Bottom row: frame 162. Vortic-ity confinement tends to create high-frequency noise everywhere, whileIVOCK produces a natural transition from laminar to turbulent flowwith realistic vortex structures along the way.to reduce the computational burden for Biot-Savart summation, in an attempt tomake these methods more practical for production. Unfortunately, these methodsare less intuitive for artists (velocity can’t be modified directly), tend to be moreexpensive (finding velocity from vorticity requires the equivalent of three Poissonsolves), and continue to pose problems in formulating good boundary conditions.We were nevertheless inspired by the mechanism to induce velocity from vorticity,and constructed our IVOCK scheme as a way to bring some of the advantages ofvortex methods to more practical velocity-pressure solvers.To balance the trade-off between pure grid-based methods and pure Lagrangianvortex methods, researchers have incorporated vorticity-derived forces in Euleriansimulations. Foremost among these is the vorticity confinement approach [54],where a force field is derived from the current vorticity field to boost it. Selleat al. [49] refined this by tracing “vortex particles”, source terms for vorticity con-finement tracked with the flow which provides better artistic control over the addedturbulence. Unfortunately, vorticity confinement relies on a non-physical parame-ter which must be tuned by the artist. As we can clearly see in figure 3.3, too smalla parameter doesn’t produce interesting motion for the early frames of a smoke ani-28mation, while still turning the smoke into incoherent noise in later frames. IVOCKon the other hand is built in a more principled way upon vortex dynamics, par-tially correcting the truncation error in velocity-pressure time-splitting, and pro-duces natural swirling motions without any parameters to tune. In addition, it isorthogonal to vorticity confinement (viewed as an art direction tool) and can hencebe used together (figure 3.13).Turbulence can also be added to the flow as a post-simulation process, such aswith wavelet turbulence [32]. In particular, the wavelet up-sampling scheme relieson a good original velocity field to produce visually pleasing animations; IVOCKis again orthogonal to this and could be adopted to enhance the basic simulation.3.3 The IVOCK schemeWhen solving Navier-Stokes, the fluid velocity is usually advanced to an interme-diate state ignoring the incompressibility constraint:u˜ = Advect(un) , (3.3)˜˜u = u˜+∆t f , (3.4)where un is the divergence-free velocity field from the previous time-step and f isa given force field which may include buoyancy, diffusion, vorticity confinement,and artistically controlled wind or motion objects.From this intermediate velocity field ˜˜u one can construct the final divergence-free velocity un+1 with pressure projection, which is usually the place where bound-ary conditions are also handled:un+1 = Proj(˜˜u). (3.5)We observe that in vortex dynamics, the intermediate velocity field u˜ of equa-tion 3.3 is analogously solved using the velocity-vorticity (u−ω) formula:1. given ωn, solveDωDt= ω ·∇u (3.6)to get ω˜;292. deduce the intermediate velocity field u˜ from ω˜ .Equation 3.6 is the advection-stretching equation for vorticity in 3D, and step2 of this (u−ω) formula is usually solved using a Biot-Savart summation, orequivalently by finding a streamfunction Ψ (vector-valued in 3D), which satisfies∇2Ψ=−ω , and readily obtaining the velocity u˜ from u˜ = ∇×Ψ .Equation 3.6 suggests the post-advection vorticity field ω∗ =∇× u˜ should ide-ally equal the stretched and advected vorticity field ω˜ , but due to the simple natureof self-advection of velocity ignoring pressure, there will be an error related to thetime step size.We therefore define a vorticity correction δω = ω˜−ω∗, from which we deducethe IVOCK velocity correction δu and add this amount to the intermediate velocityu˜. Algorithm 4 provides an outline of the IVOCK computation before we discussthe details.Algorithm 4 IVOCKAdvection(∆t,un, u˜)1: ωn← ∇×un2: ω˜ ← stretch(ωn)3: ω˜ ← advect(∆t, ω˜)4: u˜ ← advect(∆t, un)5: ω∗← ∇× u˜6: δω ← ω˜−ω∗7: δu← VelocityFromVorticity(δω)8: u˜ ← u˜+δuIn essence, IVOCK upgrades the self-advection of velocity to match self-advectionof vorticity, yet retains most of the efficiency and all of the flexibility of a velocity-pressure simulator. In particular, IVOCK doesn’t change the pressure computation(as the divergence of the curl of the correcting streamfunction is identically zero,hence the right-hand-side of the pressure projection is unchanged by IVOCK), butsimply improves the resolution of rotational motion for large time steps, which weearlier saw was hurt by velocity self-advection.30Figure 3.4: Vorticity and streamfunction components are stored in a stag-gered fashion on cell edges in 3D (red line), while velocity componentsare stored on face centers. This permits a natural curl finite differencestencil, as indicated by the orange arrows.3.3.1 Vortex dynamics on gridsData storageExtending the classic MAC grid [25] widely adopted by computer graphics re-searchers for fluid simulation, we store vorticity and streamfunction componentson cell edges in a staggered fashion, so that the curl operator can be implementedmore naturally, as illustrated in figure 3.4. For example, if the grid cell size is h,the z-component of vorticity on the z-aligned edge centred at (ih, jh,(k+ 12)h) isapproximated asωz(i, j,k+12) =1h[(v(i+0.5, j,k)− v(i−0.5, j,k))−(u(i, j+0.5,k)−u(i, j−0.5,k))]. (3.7)Vortex stretching and advection on gridIn 3D flow, when a vortex element is advected by the velocity field, it is alsostretched, changing the vorticity field.We solve equation 3.6 on an Eulerian grid using splitting; we solve∂ω∂ t= ω ·∇u, (3.8)to arrive at an intermediate vorticity field ωn+ 12 , which is then advected by the31Figure 3.5: Barnes-Hut summation [5] for the boundary values demonstratedin 2D: we construct the monopoles of tree nodes from bottom-up (left),and then for an evaluation voxel (black) in the ghost boundary region(grey), we travel down the tree, accumulating the far-field influence bythe monopoles (blue nodes), and do direct summation only with closeenough red cells (right).velocity field asDωDt= 0. (3.9)When discretizing equation 3.8 on an Eulerian grid, a geometrically-basedchoice would be to compute the update with the Jacobian matrix of u:ωn+12 = ωn+∆tJ (u)ωn. (3.10)Constructing the Jacobian matrix and computing the matrix vector multiplica-tion involves a lot of arithmetic, so we simplified this computation by usingω ·∇u = ∂u∂ω= limε→0u(x+0.5εω)−u(x−0.5εω)ε. (3.11)In our computation, with h the grid cell width, we choose ε = h‖ω‖2 , which canbe seen as sampling the velocity at the two ends of a grid-cell-sized vortex segmentand evaluating how this segment is stretched by the velocity field.Once the vorticity field is stretched, it is advected by any chosen scheme intro-duced in §3.2 to get ω˜ .32Obtaining δu from δωTo deduce the velocity correction from the vorticity difference, we solve the Pois-son equation for the stream function:∇2Ψ=−δω (3.12)We solve each component of the vector stream function separately: Let sub-script x indicate the x-component of the vector function, we have∇2Ψx =−δωx, (3.13)The equations for y and z components can be obtained in a similar fashion.In vortex dynamics, it is important for this Poisson equation to be solved inopen space to get natural motion, without artifacts from the edges of the finite grid.The open space solution of Poisson’s equation can be found by a kernel summationwith the corresponding fundamental solution. Consider the scalar Poisson problem,∇2Ψx =−δωxΨx(p) = 0, p→ ∞,(3.14)Its solution, in a discrete sense, is given by the N-body summationΨx(pi)≈ ∑all j, j 6=iδωx jv j4piri, j, (3.15)where δωx j, v j are the vortex source strength (component-wise) and volume ofeach sample point (voxel), respectively, pi is the evaluation position, and ri j is thedistance between the ith and jth sample positions. With M evaluation points (onthe six faces of the domain boundary) and N source points (the number of voxels),the direct N-body summation has an O(MN) cost. However, once an evaluationposition pi is far from a cloud of s sources, the summation of the influence of eachindividual source can be accurately approximated by the influence of the monopole33of those sources:s∑j=1δωx jv j4pi‖pi− p j‖ ≈14pi‖pi− pc‖s∑j=1δωx jv j, (3.16)where pc is the center of mass of the cloud, and ∑sj=1 δωx jv j is the so-calledmonopole. Defining an octree on the voxels, akin to multigrid, we construct frombottom-up the monopoles of clusters of voxels, and then for each evaluation po-sition, we traverse the tree top-down as far as needed to accurately and efficientlyapproximate the N-body summation, as illustrated in Fig. 3.5. Similarly to Liu andDoorly [36], we apply this summation scheme at the M ∝ N23 boundary cells toapproximate the open-space boundary values at an O(N23 log(N)) cost, which wethen use as boundary conditions on a more conventional Poisson solve.We discretize the Poisson equation with a standard 7-point finite differencestencil. Because the vorticity difference is small (one can view it as a truncation er-ror proportional to the time step), we can get away with applying a single multigridV-cycle to solve for each component of the stream function. Recall this is not thestream function for the complete velocity field of the flow, just the much smaller ve-locity correction to the intermediate self-advected velocity! One V-cycle providesadequate accuracy and global coupling across the grid at a computational cost ofonly about 20 basic red-black Gauss-Seidel iterations. This is cheaper than thepressure solve, which requires more than three V-cycles for the required accuracy,and is also substantially cheaper than a vorticity-velocity solver which requiresthree accurate Poisson solves.Once the approximate streamfunction is determined, the velocity correction iscomputed as δu = ∇×Ψ.3.3.2 DiscussionVortex stretchingBefore proceeding to the 3D applications of the IVOCK scheme for smoke §3.4.1,water §3.4.2 and combustion §3.4.3, we present a few examples illustrating theeffect of vortex stretching.34Figure 3.6: Two sequences from 3D rising smoke simulations. Top row: MC-IVOCK with vortex stretching. Bottom row: MC-IVOCK with vortexstretching switched off. With vortex stretching, vortex rings changetheir radius under the influence of other vortex rings, these process caneasily perturb the shape of vortex rings, breaking them to form newvortex structures, which brings rich turbulence into the flow field.In 2D flows, vortex stretching doesn’t take place, hence one need only advectthe vorticity field. Figure 3.7 compares SF and SF-IVOCK in 2D.In 3D flows, vortex stretching plays an important role. Rising vortex ringsleapfrogging each other is a physically unstable structure: the vortex rings breakup under small perturbations and form new vortex structures. However, due to theadded complexity of implementing the vortex stretching term, and earlier questionsabout it stability, practitioners are often tempted to simply drop it, despite this notbeing consistent with the fluid equations. Figure 3.6 illustrates that IVOCK withoutvortex stretching produces, at large time steps, visibly wrong results. In figure 3.8,we show that the ability of the IVOCK method to capture vortex stretching canbring more turbulence to the flow, such as wrinkles on the smoke front, again withrelatively large time steps.Additional forcesIVOCK is only a correction to the self-advection part of a standard velocity-pressuresolver. As such, we do not need to take into account how additional force termssuch as buoyancy, viscosity, and artist controls will interact with vorticity. Theseforces are incorporated into the velocity directly in a different step.35Figure 3.7: Left: 2D buoyancy flow simulated with SF. Right: the same withSF-IVOCK. The resolution and time step were the same; the IVOCKscheme produces a more richly detailed result.Figure 3.8: Vortex stretching enhances vortical motion, captured more accu-rately with IVOCK. Transparent renderings illustrate the internal struc-tures. Left and middle-lfet: FLIP. Middle-right and right most: FLIP-IVOCK.Accuracy vs. performanceWhile more multigrid V-cycles can be performed to improve the accuracy of thestream function solver, we observe the visual improvement is minor. In our ex-periments, more than 60% extra computation time was required for the solver toconverge to a relative residual reduction of 10−6, which is probably not worthwhile.Fig. 3.9 shows two simulations obtained with different convergence criteria.Figure 3.9: Left: Frame 150 of an IVOCK simulation with only one multigridV-cycle. Right: Frame 150 of an IVOCK simulation where multigrid V-cycles are taken reduce the residual by 10−6.36Boundary conditionsWhile solving IVOCK dynamics, we don’t take boundaries into account at all: wepose the vorticity-velocity equation in open space, and allow the subsequent pres-sure projection to handle boundaries. However, in some simulations with a strongshear layer along a viscous boundary, we found IVOCK could cause instability;this was cured easily by zeroing out δω within a few grid cells of the boundary.We hypothesize the root cause of this problem is that the inviscid fluid equationswe numerically solve predict free-slip boundary conditions at solids, allowing es-sentially an arbitrary O(1) jump in tangential velocity across a solid boundary;evaluating vorticity with a stencil which also includes a voxel with solid velocitythen results in an unphysical O(1/h) singular spike, causing instability. Unlessviscosity is taken into account and the viscous boundary layer is sufficiently re-solved on the grid, evaluating vorticity near a solid is almost certainly numericallyill-advised.It is worth noting that velocity-pressure solvers in general, and IVOCK in par-ticular, rely on the numerical diffusion of the advection to produce vortex sheddingalong boundaries (which are otherwise not predicted by the fully inviscid equa-tions). Figure 3.13 shows an example where rising smoke collides with an obstacle;the velocity boundary condition is handled by the pressure solver, and no specialtreatment is needed for the IVOCK streamfunction computation.3.4 Applications and results3.4.1 SmokeThe clearest application of IVOCK is in enriching smoke simulations, where vortexfeatures are of crucial importance visually.In 2D, we initialized a raising vortex pair from a heat source. With SF, thisvortex pair all but vanishes after 100 time steps; SF-IVOCK preserves the vortic-ity much better, as shown in figure 3.10. We recorded the total vortex strength(enstrophy) at each time step, plotted in figure 3.11.In 3D, we applied IVOCK to different advection schemes. Figure 3.1 illus-trates the qualitative improvements to all of them. For large time steps, IVOCK37Figure 3.10: Rising vortex pair initialized by a heat source. Left column: SFwith ∆t = 0.01. Middle column: SF with ∆t = 0.0025. Right column:SF-IVOCK with ∆t = 0.01. Notice IVOCK preserves vorticity bestand produces the highest final position among the three.significantly enhances the rotational motion and structures before turbulence fullydevelops.0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.800.20.40.60.811.21.41.61.8time∫ Ω||ω|| 2 dxVorticity−time curveIVOCKstab dt = 0.01Advectstab dt = 0.01Advectstab dt = 0.0025Buoyancy forceremovedFigure 3.11: Vorticity vs. time curve of the 2D vortex pair simulation. WhileIVOCK does not conserve vorticity exactly due to approximations inadvection and the streamfunction computation, it still preserves signif-icantly more.Computationally speaking, the overhead for IVOCK includes the stretchingcomputation, vorticity evaluation, vorticity advection, and three V-cycles to ap-proximate the vorticity streamfunction. These operations add about 10% to 25%extra runtime to the original method as a whole: see table 3.2 for sample per-38Figure 3.12: IVOCK can be cheaper and higher quality then taking smalltime steps. Left and mid-left: BFECC. Mid-right and right: BFECC-IVOCK with twice as large a time step, computed in much less time.formance numbers. Alternatively put, the runtime overhead of IVOCK could beequivalent to improving the original method with about 0.94×∆x in spatial resolu-tion, or taking about 0.83× smaller time steps. Figure 3.12 illustrates that runningIVOCK produces better results than even a simulation with half the time step size,despite running much faster.Table 3.2: Performance comparison of IVOCK augmenting differentschemes, for a smoke simulation at 128x256x128 grid resolution, run-ning on an Intel(R) Core(TM) i7-3630QM 2.40GHz CPU.Method sec /time step % overheadSF 26.6SF-IVOCK 30.5 14%SL3 27.7SL3-IVOCK 32.6 17%BFECC 31.8BFECC-IVOCK 39.9 28%MC 27.9MC-IVOCK 35.4 26%FLIP 45.5FLIP-IVOCK 51.4 12%3.4.2 LiquidsWe also implemented IVOCK for liquids solved as free surface flow (see figure3.14), with a solver based on Batty et al.’s [6]. In this case we only compute thevorticity correction δω inside the liquid, sufficiently below the free surface so that39Figure 3.13: Rising smoke hitting a spherical obstacle: the velocity bound-ary condition is handled entirely by the pressure solve, and doesn’tenter into the IVOCK scheme at all. This example includes a small ad-ditional amount of vorticity confinement to illustrate how the methodscan be trivially combined.extrapolated velocities are not involved in the vorticity stencil, as otherwise wefound stability issues. We solve for the velocity correction with a global solve asbefore, disregarding boundaries.Interestingly, we found IVOCK made little difference to the common waterscenarios we tested, even ones where we made sure to generate visible vortexstructures (e.g. trailing from a paddle pushing through water). We note first thattypically the only visible part of a liquid simulation is the free surface itself, un-like smoke and fire, so interesting vortex motions under the surface are relativelyunimportant. We also hypothesize that most visually interesting liquid phenomenahas to do with largely ballistic motion (splashes) or irrotional motion (as in oceanwaves), chiefly determined by gravity, so internal vorticity is of lesser importance.While potential flow can be modeled with a vortex sheet along the free surface, weleave capturing this stably in IVOCK to future work.3.4.3 FireFor combustion, the velocity field is not always divergence-free due to expansionsfrom chemical reactions and intense heating. We subsequently modify the vortexstretching term, noting that ω/ρ is the appropriate quantity to track [58]:ddt(ωρ)=(ωρ)·∇u. (3.17)40Figure 3.14: IVOCK applied to liquids. Top row: dam break simulations ob-tained with FLIP (left) and FLIP-IVOCK (right). Bottom row: movingpaddle in a water tank, simulated with FLIP (left) and FLIP-IVOCK(right). In these and several other cases we tested, FLIP-IVOCK isnot a significant improvement, presumably because interior vorticityis either not present (irrotional flow) or not visually important.We discretize equation 3.17 with Forward Euler,1∆t(ωn+1ρn+1− ωnρn)=(ωnρn)·∇un (3.18)which implies thatωn+1 =(ρn+1ρn)(ωn+∆tωn ·∇un) . (3.19)Therefore the new strength of a vortex element, after being stretched and advected,shall be scaled by ρn+1ρn .In compressible flow, mass conservation can be written as∂ρ∂ t+∇ · (ρu) = DρDt+ρ∇ ·u = 0. (3.20)Rewriting equation 3.20 gives1ρDρDt=DlnρDt=−∇ ·u, (3.21)41and by using Forward Euler again we getρn+1ρn= exp(−∆t∇ ·u) . (3.22)Plugging equation 3.22 into equation 3.17, we observe that in compressibleflow, the vorticity field after being stretched and advected, should be scaled by afactor of exp(−∆t∇ ·u). This scaling is the only change we make to IVOCK forcompressible flow. We combined the modified IVOCK scheme with traditionalvolumetric combustion models described by Feldman et al. [17]; figure 3.15 showsexample frames from such an animation.3.5 ConclusionWe argue IVOCK is an interesting stand-alone method to cheaply enrich the highlyflexible grid-based velocity-pressure framework with much better resolved vortic-ity at large time steps, and which can be applied to a variety of advection schemesand fluid phenomena. We believe it brings many of the advantages of vortex solversto velocity-pressure schemes, but with only 10∼ 25% extra computation and with-out the complexity of handling boundary conditions etc. in a vorticity formulation.However, there are some limitations and areas for future work we would highlight.Currently IVOCK is limited to fluid simulations on uniform grids. A versionfor adaptive grid or tetrahedral mesh solvers doesn’t appear intrinsically difficult,as it could reuse the solver’s Poisson infrastructure. Fully Lagrangian velocity-pressure particle methods, especially those that already require some form of vor-ticity confinement to combat intrinsic numerical smoothing, such as position-basedfluids [38], may also benefit from the IVOCK idea.The presented IVOCK scheme doesn’t exactly recover total circulation (or en-ergy) as Mullen et al.’s symplectic integrator does [42], nor does it reproduce tur-bulent flow quite as well as Lagrangian vortex methods. We also have at presentlittle more than numerical evidence to argue for its stability, and we haven’t thor-oughly investigated the trade-offs of using just an approximate multigrid solutionfor the vorticity-velocity step. It could be worthwhile to enhance the stability byexploring Elcott et al.’s circulation preserving method [14].42Figure 3.15: Fire simulations making use of our combustion model. Top row:BFECC. Bottom row: BFECC augmented with IVOCK using SL3 forvorticity and scalar fields. The temperature is visualized simply bycolouring the smoke according to temperature.The IVOCK scheme alone doesn’t provide any insight into boundary layermodels or unresolved subgrid details, hence one possible future project could be tocombine it with subgrid turbulence models, for post-processing [32] or on-the-flysynthesis [45]. However, we did observe stability problems with IVOCK compen-sation active in strongly shearing boundary layers; until a better understanding ofboundary layers is reached, we recommend only including the vorticity differencedriving IVOCK at some distance away from boundary layers.43Chapter 4A PPPM Fast SummationMethod for Fluids and BeyondFigure 4.1: Vortices of raising smoke hitting a moving ball (ball not ren-dered).4.1 IntroductionA common numerical problem in and outside of graphics is solving the Poissonequation ∇2 f = −ρ whose exact solution (assuming zero boundary conditions atinfinity for now) can be expressed as an integral with the fundamental solutionkernel, in 3D:44Table 4.1: Common symbols used throughout the paper.xi position of ith Lagrangian elementX position of a grid cell’s centeru(x) velocity evaluated for Lagrangian elementsu(X) velocity evaluated at grid cellsu f ar the velocity introduced by far field vorticityunear the velocity introduced by near field vorticityuω rotational part of the flow velocityuφ irrotational, divergence free part of the flow velocityu∞ velocity at infinity, such as the speed at which thereference frame is movingω 3D vector valued vorticity strengthσ scalar valued charge density per unit areah width of the grid cellK K = 0,1,2,3, ... the local correction rangeΨ a vector valued stream functionΦ scalar valued potential functionωL the local vorticity field defined on nearby cellsΨL the local stream function with source term ωLuL velocity introduced by ωLf (x) =ˆρ(x′)14pi‖x− x′‖dx′. (4.1)(In other dimensions, the kernel is of course a little different).In many physics applications, ρ(x) is a sum of point sources at positions x j andstrengths ρ j, or can be approximated as such: ρ(x) = ∑ j ρ jδ (x− x j) (e.g. vortexblobs [9]). In astrophysics, calculating the gravitational force by summing overall masses is essentially solving the Poisson problem with the density distributionas the right-hand side [27]. The same Poisson problem appears in electrostaticswith a charge distribution on particles (or charge panels), with a recent interestingapplication in graphics [63]. A more thorough review regarding the Poisson prob-lem and N-body dynamics can be found in Hockney and Eastwood’s work [27]. Inthose simulations, the solution is usually obtained by45f (x) = ∑j:x 6=x jρ jv j4pi‖x− x j‖ . (4.2)With N particles and M locations to evaluate, often at the N particles them-selves, this takes O(MN) time, which scales poorly.The Poisson problem or its solution as a summation over particles has appearedin diverse graphics papers. For example, Kazhdan et al. [30] solve a Poisson prob-lem to reconstruct a surface from point samples; Jacobson et al. [28] segment theinside from the outside of triangle soup by constructing a harmonic function fromsources on triangle faces, again using a summation approach; many fluid solvers(e.g. [16, 53]) solve a Poisson problem every time step for pressure projection; andLagrangian vorticity-based fluid solvers (e.g. [44]) evaluate the velocity field bysumming over vortex sources with the Biot-Savart law, which is just the curl ofequation (4.2).While sometimes it is natural to discretize the Poisson equation on a grid withfinite differences, and impose some boundary condition on the edges of the grid,often we want to solve the “free space” problem with boundaries only at infinity,or with sources or boundaries distributed arbitrarily through space: on particles, oncontours, on the triangles of a mesh, etc. Here the integral or summation approachwould shine but for the cost of directly evaluating the sum.Alternative methods exist to accelerate the expensive summation. Foremostamong these is the Fast Multipole Method (FMM) [23], which can reduce the costfrom O(MN) to O(M+N). Probably due to implementation complexity and (typi-cally) a large overhead in runtime, the FMM has not been widely adopted in graph-ics: parallelized direct summation can easily outperform FMM for N up to 100,000in our experience.Vortex-in-cell (VIC) methods [12] have also been adopted to efficiently al-beit roughly approximate the velocity field. VIC splats the right-hand side of thePoisson problem to a fine grid, solves the PDE on the grid with finite differencesor similar, and interpolates the solution back from the grid. The computationalcomplexity is typically dominated by an efficient grid solve. However, VIC losesaccuracy when the vorticity fields are under-resolved on the grid resolution; trans-ferring values back and forth between grid and particles forms a major source of46numerical diffusion.A more promising approach, in our view, is the Particle-Particle Particle-Mesh(PPPM) method [1]. The Particle-Mesh component of PPPM approximates thesmooth long-range interactions by splatting the particles on to a grid (mesh), usinga fast solver for a finite-difference version of the Poisson problem on that grid, andthen interpolating back to evaluation points, just like VIC. The Particle-Particlecomponent greatly enhances the accuracy by including the kernel summation justfor very close pairs of particles, whose interactions aren’t properly captured withthe grid. Assuming the distribution of particles is such that a moderate resolu-tion grid exists with a small number of particles per cell, the particle operationsin PPPM are O(M+N) and the overall runtime is determined by the grid solver(which traditionally has been O(N logN) with an FFT solver). Moreover, opera-tions such as splatting to or interpolating from a grid are easy to implement withlittle overhead. That said, prior PPPM methods have had trouble achieving linearruntime (due to the fast solver used) and full accuracy (due to boundary condi-tions at the edges of the grid, and errors associated with near-field interactions ascomputed on the grid).This paper first extends PPPM with two significant algorithm improvements:• a new “monopole” grid boundary condition to accurately simulate unboundeddomains,• and a new error correction for nearby particle contributions in the grid solvethat is simple, fast, and accurate.We also upgrade the traditional fast solver to linear-time multigrid, and parallelizeit all on the GPU for exceptional performance.In the second part we turn to applications of our new PPPM, in particular anew vortex particle smoke solver. Lagrangian vortex methods succinctly representrichly detailed flow, and suffer no numerical diffusion or projection-related dis-sipation as other solvers may. Their central step is reconstructing velocity fromvorticity via the Biot-Savart law, essentially the curl of equation (4.2),u(x) = ∑j:x 6=x j∇14pi‖x− x j‖ ×ω j (4.3)47or equivalently solving a Poisson problem for a potential Ψ and taking its curlto get the velocity. Inviscid solid boundary conditions can be implemented via apotental flow correction to this velocity field, which through a boundary integraldiscretization with an iterative linear solve likewise relies on an N-body summation/ Poisson solve. We use our new PPPM method in each case to achieve a highperformance, high quality smoke simulation. The greater accuracy of our PPPMapproach gives visually superior results than simpler techniques such as Vortex-in-Cell and truncated kernels.Apart from the challenge of efficient N-body summation, prior vortex particlemethods in graphics have had difficulties with the vortex stretching term (necessaryfor conservation of angular momentum and the characteristic look of developingturbulence) and vortex shedding from solid boundaries. We propose solutions:• an easy-to-implement vortex-segment-derived treatment of vortex stretchingwith automatic conservation properties,• and a simple boundary layer treatment to shed vortices from solids, basedon injecting vortex particles to boost the potential flow no-stick boundarycondition to a no-slip boundary.We also show that PPPM can be exploited for non-physics-related problems ingraphics, such as Poisson surface reconstruction.4.2 Related workVortex methods are popular for high Reynolds number computational fluid dynam-ics [11], very efficiently capturing turbulent details; as vorticity in particular is vi-sually very important, many graphics researchers have adopted vortex approaches.Yaeger et al. [65] began the trend with a VIC solver to produce the realistic lookof atmospheric flow on Jupiter. Angelidis and Neyret [4] used Lagrangian vortexprimitives for their flexibility; Selle et al. [49] used vortex particles to augment amore traditional Eulerian velocity-pressure simulation. Pfaff et al. [45] used vortexparticles from pre-processed boundary layers to synthesize turbulence. Weißmannand Pinkall [64] used Biot-Savart summation for vortex rings.48Brochu et al. [8] captured vorticity with a Lagrangian mesh, including its gener-ation from buoyancy, and used FMM for linear-time Biot-Savart summation. Pfaffet al. [46] used the same vortex sheet representation to generate and evolve smallscale details, but accelerated summation with a truncated kernel, relying on a coarsegrid simulation to account for global effects. Pfaff et al.’s hybrid use of a grid andnearfield vortex summation is reminiscent of PPPM, but lacks the error correctionstage we develop, and suffers from the numerical dissipation of the pressure pro-jection on the grid. Both Brochu et al. and Pfaff et al. also pay a high cost in meshoperations to maintain the deforming vortex sheet, which we avoid by trackingdensity and calculating buoyancy with unconnected particles.Despite their attraction for turbulence, vortex methods have not been widelyaopted in the graphics industry, perhaps because of the difficulty in acceleratingBiot-Savart summation, or boundary condition challenges, or the intricacies of vor-tex shedding. Addressing these problem motivates this paper.Many graphics applications outside fluids also solve the Poisson problem. InPoisson surface reconstruction [30], the solution is important only around a narrowband of the input, but the associated Poisson problem is most naturally posed onan infinite domain, making an awkward fit for traditional voxel approaches.Gumerov and Duraiswami [24] demonstrated a high-performance GPU imple-ment of FMM, reporting 1.4s to sum N = 106 particles on an NVIDIA GeForce8800 GTX. Our PPPM code runs the same problem on an NVIDIA GeForce GT650M (laptop) in 0.7s. More subjectively, our experience is that implementing afast PPPM is vastly easier than even basic FMM.We use multigrid as our fast grid solver as have many other graphics papers(e.g. [18, 40, 41]), though we should note that Henderson suggests on shared mem-ory multiprocessors FFT still may be faster [26].4.3 Particle-particle particle-mesh method(PPPM)In this section we take Biot-Savart summation as an example to explain the PPPMalgorithm in detail. See figure 4.2 for an overview.Given the vorticityω =∇×u of an incompressible flow, and ignoring boundaryconditions for now, we can reconstruct velocity (up to addition of the gradient of49Find bounding box of vortexi i f rt Determine the global domaint r i t l l iAssign particle values to grid, compute B.C.i rti l l t ri , t . .Solve Poisson system to get global velocity fieldl i t t t l l l it fi lFor each gridcell, cancel the short range contributionr ri ll, l t rt r tri tiFar field constructionuglobalufar ufarufinalVelocity evaluationGet far field velocityby interpolationt f r fi l l it i t r l tiAdd near field direct sum-mation to get final velocity r fi l ir t -ti t t fi l l itFigure 4.2: An overview of our PPPM algorithm, which consists of a far-fieldconstruction step and a velocity evaluation step.a harmonic scalar potential) by solving for a streamfunction Ψ whose curl is thevelocity: ∇2Ψ=−ωΨ(x) = 0 x→ ∞ (4.4)In the case where the right hand side of this Poisson system is given by a sumof n vortex elements, each carrying a vortex strength of ω i, Ψ can be found bysummation (c.f. equation 4.2). The velocity is just the curl of this sum (equation4.3), further expressed asu(xi) =14pin∑j=1, j 6=iω j× (xi− x j)∣∣xi− x j∣∣3 , (4.5)which is known as Biot-Savart summation.When evaluating at a specific point, we decompose this velocity as u = u f ar +unear, where unear is the contribution from nearby vortices and u f ar gathers influ-ences from the rest.Due to the singularity of the kernel, unear varies rapidly in space and is relatedto the small scale motion (turbulence), while u f ar is smooth and captures the largescale motion.Realizing that multigrid Poisson solvers are excellent for smooth, large scalemotions, while particle-particle direct summations are promising for the turbulentsmall scale motions, PPPM first estimates u f ar for each grid cell using a particle-mesh solver (§4.3.1) and local-cancellation (§4.3.2). After that, PPPM interpolates50u f ar for each vortex element and evaluates unear using direct summation. Detailsof the computational steps in the PPPM Biot-Savart summation are outlined inAlgorithm 5, where the subscript p indicates the quantity is carried by a particle,and the subscript g indicates the quantity is defined on a grid.Algorithm 5 PPPMBiotSavartSummation(ω p, vp, xp, N, up)1: BB←GetBoundingBox(xp)2: DetermineComputationDomain //3×BB3: ωg←ParticleToMesh(ω p)4: ComputeMonopoleDirichletBoundaryCondition5: Ψ←ParallelMultigridSolvePoisson6: ug←Curl(Ψ)7: for each grid cell in parallel8: Subtract near field estimate to get far field component9: endfor10: for each particle i in parallel11: up,i← InterpolateFarField12: endfor13: for each particle i in parallel14: Sum over particles in neighborhood η(i) of i:(evaluating an accurate near field component)15: ∆u← ∑ j∈η(i) ω p, j×(xp,i−xp, j)4pi|xp,i−xp, j|3(1− exp( |x−x j|α)3)16: up,i← up,i+∆u17: endfor4.3.1 Particle-mesh step in open spaceIn the particle-mesh step, we determine the large scale flow motion by solving aPoisson equation with Dirichlet condition(on domain boundary):∇2Ψ=−ω in ΩΨ= g on ∂Ω (4.6)The imposition of an artificial boundary at the edge of the grid is necessaryfor the solve, but makes approximation of the open space (infinite domain) prob-lem tricky. We propose two special advances in our discretization to gain higher51accuracy for the large scale motion in open space, compared to prior PPPM meth-ods which use periodic boundaries or a homogenous Dirichlet condition. First, weuse a computational domain that is three times as large as the bounding box of thevortex particles, giving us an effective padding region: the artificial boundary isdistant from all sources in the problem. In open space the true solution Ψ at anyevaluation point xb on the boundary is exactlyΨ(xb) =N∑i=1ω (xi)h34pi |xi− xb| . (4.7)Computing these values to be used as the Dirichlet boundary condition g wouldgive, up to truncation error, the exact open space solution. However evaluatingequation 4.7 directly is too costly: instead, we compute the cheap monopole ap-proximation of the particle quantities, replacing equation 4.7 withΨ(xb) =N∑i=1ω (xi)h34pi |xi− xb| ≈∑Ni=1ω (xi)h34pi |xc− xb| =m(xc)4pi |xc− xb| (4.8)where xc is the center of all the vortex particles. The monopole is accurate whenall the vortex particles are far away from the domain boundary, which we guaranteeby our domain construction. We dub this the “monopole boundary condition” (line4 in Algorithm 5), which experimentally we found doubles the accuracy of PPPM.After we set up the computation domain and construct the boundary condition,the particle values are splatted to the grid usingω (X) =1h3∑iviω (xi)(∏θ=1,2,3wθ(Xθ − xi,θh))(4.9)where subscript θ = 1,2,3 indicates the corresponding component of a 3Dvector, h the size of the grid cell, v the volume of a vortex blob, i the index of avortex blob, and w the splat kernel.For the choice of w, we simply use nearest grid point (NGP) assignment, mean-ing each grid cell gathers the total strength of the vortex particles inside it:wθ (r) =1, r ∈ [−0.5,0.5) ,0, else, (4.10)52In our implementation, we accelerate the process by parallelizing for each gridcell to gather vortex values around it, using spatial hashing [59] for efficient neigh-bour finding, following Green’s implementation [22].We then discretize the vector Poisson system using the seven-point secondorder finite difference scheme, arriving at three scalar Poisson systems, one foreach component. We solve the resulting linear system for the streamfunction usingmultigrid, following Cohen et al.’s implementation [10].Once we have found the streamfunction, we take its curl to get a velocity fieldupm. upm is then used to derive the far field approximation for each grid cell, usingthe procedure in the next section.4.3.2 Cancelling local influences in the gridTo cancel the (relatively inaccurate) local contributions from nearby grid cells, weneed to estimate a local stream function ΨL that solves∇2ΨL =−ωL (4.11)in open space, where ωL is the near-field vorticity only from the nearby gridcells.Conceptually, we can achieve this by puttingωL back in the global domain(grid),and solving the Poisson system again. (However, as we need a different estimatefor each grid cell containing particles, in practice this would be far too expensive.)More formally, putting ωL back to the global domain is a prolongation,PωL.The global Poisson solve can be denoted as Ψ =L −1PωL. Reading back ΨL isa restriction of Ψ, R(L −1PωL). Finally, the local velocity field due to ωL isreadilyR(L −1PωL).The accuracy and efficiency of local cancellation depends on fast localizedsolutions for ΨL, say, expressed as a linear operator A:ΨL = Ah2ωL, (4.12)with A≈ 1h2RL −1P .Theuns [60] approximated A with an open space Green’s function. however,53this is quite different from the inverse of the grid-based operator, as the Green’sfunction is singular at r→ 0 where the grid-based inverse is bounded; this rendersit unable to evaluate the local contribution made by a grid value to itself. Walther[62] proposed a more accurate estimate at the cost of a pre-computation stage: witha local correction window of size K, the method stores a O(K6)matrix A whichsolves ΨL = AωL. This computation is coupled to the grid resolution, restrictingthe simulation to static grids.We instead make entries of A dimensionless quantities, allowing us to adaptthe computation domain with more flexibility. We uncover a more accurate rela-tionship between A and the open space Green’s function. For instance, with theGreen’s function approach, one hasΨL =Gh2ωL (4.13)where Gi, j = h4pi|Xi−X j| for i 6= j. (Observe that∣∣Xi−X j∣∣ is of order h, makingGi, j dimensionless.)Our observation is that the off-diagonal coefficients Ai, j are very close to Gi, j,the open space Green’s function evaluated at different grid cell centres. For diago-nal terms, instead of being 0, we find Ai,i→ 0.25 for large domains. This diagonalconstant responds to the contribution made by a grid cell to itself.Furthermore, if we replace the diagonal of G with this constant and computeits inverse, the seven-point finite difference stencil is essentially revealed. Whilewe do not yet have a full derivation, we provide our evidence in figure 4.3. Thesenumerical findings give us a new formulation to compute ΨL from ωL. We outlinethis procedure in Algorithm 6.After we obtain ΨL, we can proceed to compute uL by taking the finite differ-ence curl of ΨLuL = ∇h×ΨL (4.14)Following this naively requires each grid cell to have a small buffer for thespatial varying function ΨL, making parallelizing impractical for GPU’s memorylimitation. Instead, we return to the construction ofΨL, first looking at off-diagonalcontributions:540204060801001201400 20 40 60 80 100 120 1400.050.10.150.20.25Visualization of −AL−10204060801001201400 20 40 60 80 100 120 1400.050.10.150.20.25visualization of G0 20 40 60 80 100 120 140−6−5−4−3−2−101Row 64 of G−1Figure 4.3: Relationship between inverse of finite difference operator andGreen’s function. Left: the local inverse matrix defined by A. Mid-dle: the local inverse matrix G constructed using Green’s function andthe diagonal terms from A. Right: one row of the inverse of G, almostrevealed the 7-point stencil.Algorithm 6 Local_inverse(ωL)1: ΨL← 02: for i, j,k ∈ L //L = {i2, j2,k2|i2 ∈ [i−K, i+K], j2 ∈ [ j−K, j+K],k2 ∈ [k−K,k+K]}3: X1← X i, j,k4: if i2, j2,k2 ∈ L, i2, j2,k2 6= i, j,k5: X2← X i2, j2,k26: ΨL (X1)←ΨL (X1)+ h3ωL(X2)4pi|X1−X2|//using Green’s function7:8: else //i2, j2,k2 = i, j,k9: ΨL (X1)←ΨL (X1)+0.25h2ωL (X2)//notice X2 = X1 where, 0.25h2 is thediagonal constant10: endelse11: endfor∇X ,h×ΨL = ∇X ,h× ∑X ′ 6=Xh3ωL (X ′)4pi |X−X ′|= ∑X 6=X ′h3ωL(X ′)×(−∇X ,h 14pi |X−X ′|)(4.15)where ∇X ,h 14pi|X−X j| denotes a discrete gradient operator to the Green’s func-55tion at X . Let e1 be the unit vector (1,0,0); we have(∂∂x)X ,h14pi |X−X ′|=12h(14pi |X +he1−X ′| −14pi |X−he1−X ′|). (4.16)When X±hek−X ′ = 0, 14pi|X+hek−X ′| is set to be 0.We then compute a streamfunction buffer from the diagonal contributions,Ψd = 0.25h2ωg (4.17)We then take the discrete curl of Ψd to get ud . The velocity introduced by thenear field grid is consequently obtained viauL (X)=∑j 6=ih3ωL(X ′)×(−∇X ,h 14pi |X−X ′|)+ud (4.18)The new procedure overcomes the memory overflow issues with the naive ap-proach, is a lot faster and more scalable, while giving exact same output (up toround-off) as the direct implementation.4.3.3 Velocity evaluationThe far field influence at any given grid position X can be readily obtained via:u f ar (X) = upm (X)−uL (X) (4.19)where uL is determined by equation 4.18 in §4.3.2.To evaluate the velocity for an arbitrary evaluation position xi, we first interpo-late from the far field buffer using trilinear interpolation,u f ar (x)← TriInterpolate(u f ar, x) (4.20)56and then select all the nearby vortex elements within a cube CL of size (2K+1)hcentred at x for near-field direct summation.u(x) = u f ar+ ∑x j∈CL,x j 6=xω j× (x− x j)4pi∣∣x− x j∣∣31− exp(∣∣x− x j∣∣α)3 (4.21)4.3.4 PPPM discussionFor analysis, we assume a reasonably uniform distribution of vortex particles inthe computational domain. In the case where vortices were initialized on a surfacesheet, the turbulent motion quickly breaks this sheet into many small blobs to makethe vorticity distribution roughly uniform.We divide the domain into equal sized grid cells so that each cell contains sparticles on average, hence the number of cells is proportional to Ns . Transformingthe particle quantities to the grid takes O (N) time; applying the monopole bound-ary condition takes O(N+(Ns)2/3) operations; the multigrid Poisson solver withNs unknowns takes O(Ns)operations to get a solution; computing the discrete curlor gradient of the field requires also O(Ns)operations; interpolating from the gridquantity back to particles takes O (N) operations; and finally direct summationwith nearby particles takes O (N) time.To evaluate the velocity at some other m points, the runtime analysis is similar,only replacing the interpolating and local correction procedures in the aforemen-tioned pipeline, hence we end up with O (m+n) complexity as claimed.We tested the PPPM summation code for a random distribution of vortex parti-cles, comparing the PPPM results with direct summation, and looked at the error ina weighted average manner (particles with a vanishing velocity get small weights),As we can see from figure 4.5, the PPPM solution gets closer to direct summa-tion when larger local correction (LC) windows are used. With a window size ofK = 3, we obtained an error under 1%. Furthermore, since we bound the numberof particles per cell, with 16384 vortex particles the cell size h is half the size used57102 103 104 105 106 10710−410−310−210−1100101102103104number of N−bodiestimgingPerformance of the PPPM fast summation method measured performanceO(N)O(N 2)Figure 4.4: Performance of the PPPM fast summation. Computation timegrows linearly with the number of computational elements.for 2048 vortex particles. Therefore for the same K we have smaller physical LCradius in the 16384 case, but the accuracy doesn’t decrease because one has bettergrid resolution to resolve the near field physics.The performance of the PPPM summation is shown in figure 4.4. The compu-tation time in our experiment grows linearly with respect to the number of vortexparticles.In table 4.2 we use direct summation as the reference solution to measure theaccuracy of different approximations. We compare between our Monopole B.C.(MBC) and prior work’s homogeneous B.C. We can see clearly that the MBC580 0.5 1 1.5 2 2.5 310−310−210−1100Local correction range Kaverage error compare to direct summationaccuracy of PPPM summation 256 vortex particle2048 vortex particle16384 vortex particleless than 1%Figure 4.5: Accuracy statistics of the PPPM fast summation.enhances the approximation in either case, while our PPPM method, being ableto cancel the grid cell self-influence, gives approximation an order of magnitudehigher than Theuns [60]. On the other hand, the particle-mesh method alone, evenat higher cost, gives a poor approximation to direct summation.59Table 4.2: Accuracy of different method with and without the monopole B.C.Method K w/ MBC w/o MBCPPPM 643 3 0.46% 1.13%PM 1283 N/A 6.19% 6.38%PPPM 643 w/o diagonalcancellation [60] 3 2.78% 2.86%4.4 Vortex-Particle SmokeIn vortex methods, given the vorticity distribution and boundary geometry, the ve-locity at any position x can be found as [11]u(x) = uω +uφ +u∞ (4.22)which is a combination of a rotational flow uω from the vorticity, a potential flowuφ determined by the solid objects §4.4.2, and a prescribed velocity field u∞ definedby an artist (it is better for u∞ to be divergence-free, otherwise, we suggest using∇×u∞ as a force field).Ignoring viscosity, the dynamics of vorticity can be modeled from the curl ofthe inviscid momentum equation,DωDt= (ω ·∇)u+β∇ρ×g. (4.23)A vortex blob moves according to the combined velocity u of equation. 4.22, withvortex stretching((ω ·∇)u) due to flow deformation in 3D and baroclinic vorticitygeneration β∇ρ×g concentrated on the density interface.In the presence of solid objects, an irrotational divergence-free flow field uφcan be found to cancel flow penetration at solid boundaries. Apart from modellingsolid objects as boundary conditions in the vorticity solver, solid objects also canserve as a single layer source of “charge distribution” that introduces a harmonicpotential whose gradient removes boundary penetration. More details about howto determine and use this \"electrostatic\" field can be found in §4.4.2.We discretize the flow using a set of vortex blobs, with position xi, volume vi,density ρi, velocity ui, and vortex density ω i. At each time step, we compute uω60with Biot-Savart,uω (x)= ∑all j, j 6=iv jω j× (x− x j)4pi∣∣x− x j∣∣31− exp(∣∣x− x j∣∣α)3 , (4.24)using PPPM fast summation. Notice we mollify the Green’s function kernel todesingularize the summation, a common necessity for vortex particle codes.Subdividing the surface of solid objects into many small panels, the potentialpart of velocity, uφ , can be determined byuφ (x)= ∑all j, j 6=iA jσ j (x j− x)4pi∣∣x− x j∣∣31− exp(∣∣x− x j∣∣α)3 (4.25)where A j is the area of jth panel, σ j is the potential source strength per unit areafound on the single layer, and x j is the centroid of the jth panel.At each time step, the strengths and positions of vortex blobs are updated by:1. Initialize a vortex segment for each vortex blob based on the input vortexstrength. §4.4.12. Update the position of each end of the vortex segment with the velocity de-termined by equation 4.22, hence, the vortex segment is stretched automati-cally. §4.4.13. Reduce instabilities introduced by vortex stretching.§4.4.14. Add baroclinic vorticity generation to the vorticity field.5. Add vorticity due to vortex shedding.4.4.24.4.1 Vortex segment, vortex stretching and stability issueThis section describes the way to model vortex stretching using vortex segments.Vortex stretching introduces instability in some situations, which we address with61Vortex blob to vortex segment ortex blob to vortex seg ent Vortex segment being stretched ortex seg ent being stretched Back to vortex blob representationack to vortex blob representationFigure 4.6: We switch to a vortex segment representation of vortex particlesat the beginning of each time step, move both ends of the vortex segmentin the flow, then switch back to vortex blobs.a geometry-inspired filtering approach.Vortex segmentsVortex stretching is the rate-of-change of vorticity due to the deformation of fluid.This step is naturally handled by vortex ring or sheet representations: the stretchingof the geometric elements automatically captures the vortex stretching term. Thisdoes not extent to vortex particles with no geometric extent. An obvious solutioncould be updating the vortex strength usingωn+1 = ωn+ωn ·∇un (4.26)It is unwise to take this approach because evaluating the gradient of velocity re-quires the second derivative of the Green’s function. The singularity of this func-tion introduces large numerical instabilities. Furthermore, computing the dot prod-uct is expensive.We instead use a vortex-segment approach. A vortex segment is a small spin-ning cylinder whose central axis is aligned with the vorticity direction, with twoends xa and xb, length h, and constant circulation κ . A vortex blob with vorticitystrength viω i can be translated to a vortex segment of length hi, with unit directiontˆi = ω i√‖ω i‖2 , and circulation κi = ω i/(hitˆi). We translate our vortex blob into such62Vortex segment being stretched in a non-smooth way generates numerical errori ii lImprove stability by forcing the change made by stretching to be smooth without actually diffuse ω Figure 4.7: Sudden discontinuous motion of vortex segments introduces andamplifies numerical error, which is reduced then by smoothing thestretching terms.vortex segments at the beginning of each time step, evaluate the velocity accordingto equation 4.22 for each end of the vortex segment, and move each end of thevortex segment according toxn+1i,a = xni,a+∆tuni,axn+1i,b = xni,b+∆tuni,b. (4.27)When both ends of the vortex segment are updated, the vortex segment is stretchedand we transform it back to a vortex blob with vortex strength,viω i = κi(xn+1i,b − xn+1i,a). (4.28)Notice the circulation is conserved. This whole process is illustrated in figure 4.6.Stability issuesVortex stretching is potentially unstable. Without addressing this sensitivity, oursimulation quickly diverges even with small time-steps. We noticed that in tur-bulent flow, even when the magnitude of the velocity is small, the gradient of thevelocity field can be very large, which makes ω ·∇u large and drives the numer-ical instabilities. To address this problem, we realized that in the physical world,the rate of change of vorticity due to vortex stretching has some smoothness. We63therefore compute the rate of change of vorticity due to stretching by(ω ·∇u)ni =κi(xn+1i,b − xn+1i,a)−ωn∆t(4.29)and then apply a Gaussian filter to smooth (ω ·∇u)n in the domain to get ˜(ω ·∇u)nifor each particle. Vorticity is updated byωn+1i = ωni +∆t ˜(ω ·∇u)ni . (4.30)This approach is effectively a geometric repairing that forces the deformation ofa vortex segment to be consistent with nearby vortex segments, as illustrated infigure 4.7. With this approach, we stabilized the simulation and reliably achievedlong simulations when using a constant large time step. Notice that in this scheme,we smooth the update instead of the quantity itself to preserve sharp features asmuch as possible, similar in spirit to FLIP [68].4.4.2 Solid boundaries and vortex sheddingWe incorporate boundary conditions by introducing an irrotational, divergence freeflow field uφ that cancels the normal velocity flux on the boundary made by uω +u∞. We define uφ as the gradient of a scalar potential function Φ, solved withboundary integral equations. After enforcing the no-through boundary condition,we compute a vortex sheet on the surface heuristically [9], and merge this surfacevorticity into the vorticity field.No-through boundary conditionThe no-through boundary condition can be enforced by solving for a scalar poten-tial field that satisfies Laplace’s equation with Neumann boundary condition[8]∆Φ = 0∂Φ∂n=(usolid−u f luid) ·n (4.31)We write this function as a single layer potential, with a continuous source64distribution σ on the solid boundary,Φ(x) =ˆ∂Ωσ (y)4pi |x− y|dS (y) . (4.32)Taking the normal derivative of Φ and substituting it into equation 4.31 gives aFredholm equation of the second kind,b(x) =−σ (x)2+ˆ∂Ωσ (y)∂∂nx14pi |x− y|dS (y) (4.33)where f (x) =(usolid (x)−u f luid (x)) ·n(x).We discretize this equation on a set of M boundary elements using mid-pointrule to arrive atbi =−σi2 +∑all j∂∂niσ jA j4pi∣∣xi− x j∣∣ (4.34)where A j is the area of j’th surface element, and xi and x j are the mid-points ofcorresponding boundary elements.In practise, this equation gives a linear system for σ that is very well-conditioned:iterative solvers like BiCGSTAB converge inO (1) iterations. However, the M×Mcoefficient matrix is dense and evaluating the matrix vector product is M2. To over-come this challenge, we reformulate equation 4.34 using the PPPM framework.Given a source distribution σ , with proper reordering, the off-diagonal part sum-mation is exactlybo f f−diagi = ni ·(∑all j 6=iσ jA j (x j− xi)4pi∣∣xi− x j∣∣3)(4.35)which takes the same form of evaluating a gravitational force based on mass parti-cles, where, the “mass” of the particle here are defined as σ jA j. Hence our PPPM-accelerated evaluation of the matrix vector multiplication directly follows the rou-tines described in Algorithm 7.PPPM is used in two different places here. During the iteration, we use PPPMto compute the force based on the current estimate (line 7 of algorithm 7), andwhen the iteration is terminated we use PPPM to compute uφ based on the single65Vorticity fomation on the surfaceorticity fo ation on the surfaceεγv tangentFigure 4.8: Vortex shedding process. In our approach, the vorticity strengthis explicitly determined.layer density.Algorithm 7 PPPM_accelerated_Ax(b, σ , p, A, n, M)1: f ← 02: b← 03: m← 04: for i=1:M in parallel5: mi = Ai ∗σi6: endfor7: f ←PPPM_Compute_Gravity(m, p, M)8: for i=1:M in parallel9: bi← ni · f i−0.5∗σi10: endforVortex sheddingThe inviscid assumption breaks down near boundaries, because it never introducesvorticity into the flow field. In reality fluid viscosity, no matter how small, gener-ates vorticity concentrated along a thin boundary layer along solid surfaces. In highReynolds number flows, this thin viscous boundary layer doesn’t affect the validityof inviscid approximation being made elsewhere in the flow, but rather serves as asource emitting vorticity into the flow.Chorin [9] modelled this process heuristically. in 2D. He assumed a constantvorticity strength on the boundary element and determined the vortex strengthbased on the tangential velocity slip. Extending this idea to 3D, we determine the66Figure 4.9: Rising smoke using different number of vortex particles. Left,2049 particles: middle, 16384 particles; right, 130K particles.surface vorticity direction based on the tangential velocity and the surface normal:vorticity is required to be perpendicular to these two vectors, its strength is deter-mined so that at a position in the normal direction, ε away from the surface, thevelocity introduced by this vorticity matches the tangential slip. In other words,if we put a vortex of strength A jγ j at the position ε away from the surface, thevelocity it introduces cancels the tangential slip at the boundary. This process isillustrated in figure 4.8.Once the concentrated vorticity strength is determined, we release it to theflow field by adding some amount of its value to the nearest vortex blob around thesurface or release them to a position at a small distance away from j’th boundaryelement. The amount to be released is determined as∆ω = ∆tcA jγ j (4.36)4.5 Results and conclusions4.5.1 PPPM for vortex flowWe implemented the PPPM algorithm on the GPU (nVidia GeForce GT 650m).As we can see in figure 4.4, the performance of PPPM scales linearly with thecomputational elements involved. On our machine, simulations are still runninginteractively with even 16K vortex particles on a single laptop GPU. We generatea preview simulation with a small amount of tracer particles (16K) at interactive67Figure 4.10: Without the far-field motion, direct summation using a cut-offkernel results in wrong animation. Left: simulation uses cut-off di-rect summation after 180 time steps. Right: simulation uses PPPMsummation after 180 time steps.Figure 4.11: Comparison of VIC and PPPM. Top row, sequence of 643 VICsimulation; bottom row: PPPM using the same resolution grid. Noticethat the large scale motion of the two simulation matches before thepoint where turbulent motion becomes dominant.rates, and then produce an enhanced result by just using more particle tracers (6M).Figure 4.9 illustrates rich turbulent phenomena even with a small number ofvortex particles.To achieve our final results, we are not simply interpolating velocity to thetracer particles, but rather computing the velocity of each tracer particles with thefull-influence Biot-Savart summation from vortex particles. It is only with fast68summation that it is feasible to produce the results in figure 4.1, where 130K vortexparticle and 6 million tracers were used. Each time step takes 100 sec to processon a laptop with nVidia’s Geforce GT 650M graphics card.We also observed that in our computation, local direct summation dominatesthe computation time. However, with only this near field direct summation (trun-cating the kernel to finite support), one obtains unrealistic animations. Figure 4.10,a cut-off direct summation uses the same cut-off range as the PPPM summationuses, takes almost same computation time as PPPM takes, but the smoke fails torise.On the other hand, the PPPM code without local correction is reduced to a VICsolver, which fails to produce small scale motions because of numerical smoothing.VIC running at higher resolution, to produce turbulent animations similar to PPPM,takes 64 times the memory cost and 10 times the simulation time every time step.Figure 4.11 shows representation frames.Counter-intuitively solving for the tracer particle motion takes more time in a2563 VIC simulation, even though this is just interpolation, whereas in 643 PPPMthere is a far-field interpolation followed by a more costly near field direct summa-tion. We suspect this is because memory access with larger velocity buffer is lessefficient.Our no-stick boundary condition and vortex shedding model handles differentboundary geometry robustly and produces visually plausible animations. We leftthe shedding coefficient as an artist controlled parameter. In the example shownin figure 4.12, we used a shedding coefficient c = 0.6, set the size of the movingobject to one unit length, and let it move at a speed of 4 unit lengths/sec. The vortexshedding model is able to produce complex turbulent wake patterns.4.5.2 Applying PPPM to Poisson surface reconstructionGiven a set of n sample points at position {x |x = xi, i = 1,2,3 . . .n}, with normalsnˆi at xi, the Poisson surface reconstruction [30] algorithm reconstructs the surfaceof the geometry in two steps:1 Seek a scalar field φ whose gradient best matches the vector field V defined69Figure 4.12: Moving objects in slightly viscous flow generate turbulentwakes. Top row: vortex shedding from a moving ball. Bottom row:a moving bunny.by the input, i.e., find φ that solves∇2φ = ∇ ·V (4.37)2 Determine a constant c s.t. the isosurface defined by φ (x) = c is a goodmatch of the geometry to be reconstructed. Here, c = 1n ∑i φ (xi) (average ofφ at input positions).A detailed discussion of this algorithm is beyond the scope of this paper. Here weemphasize the computation. In the original paper, an adaptive Poisson solve onan octree was used. In our approach, we only need to define the right-hand-sideon a narrow band of voxels near the input point clouds, and we can evaluate φ bysummation.More precisely, we obtain ∇ ·V in the narrow band, then solve for φ withφ (xi) =n∑j=1, j 6=ih3 f j4pi‖xi− x j‖ (4.38)here, f j = −(∇ ·V ) j on voxel j, xi and x j the position of i’th and j’th voxel, re-spectively. Those voxels and f j’s are then viewed as particles with mass, allowingus to use PPPM to calculate the summation.We tested PPPM surface reconstruction on a bunny and a dragon, shown infigure 4.13. While quality surface reconstruction is not the focus of this paper, weargue this shows the potential of PPPM for other branches of graphics. We areneither using super high resolution sparse voxels nor taking any effort in choosing70Figure 4.13: PPPM Poisson surface reconstruction of: left) a bunny, right) adragon.a good Gaussian kernel when splatting normals to construct the vector field V . Forthe dragon case we reconstructed, the voxel size is determined to be 1256 of thelongest dimension. Since the computation involves only those sparse voxels, thesummation approach became more suitable; it would be difficult for typical finitedifference approaches to impose useful boundary conditions on the boundary cellsof the narrow band.4.6 LimitationsUnlike FMM, where one can get a desired accuracy by taking enough multipoleexpansions, in PPPM, further accuracy can not be achieved for local correctionrange greater than 7 grid cells. Not only is the local cancellation imperfect, butalso, interpolating the far field has limits to its accuracy.The proposed PPPM focuses only on the N-body problem with particles, orboundary elements that can be viewed as particles. Extending the PPPM sum-mation framework to non-particles such as higher-order surface sheet or boundaryelements should nonetheless be straightforward. One could switch to higher-orderquadrature rules (or even exact integrals) for near field elements, or for each ele-ment generate samples based on quadrature rule, interpolate the value and scale itwith the quadrature weights.In FMM, one tracks adaptively distributed computational elements with adap-tive data structure, whereas PPPM use uniform background grids and uniformspace hash. This greatly simplifies the implementation and reduces computationaloverhead, but limits PPPM scalability for sparse data. Adaptive Poisson solvers[37] might address this problem.714.7 Future workAlgorithmically, many improvements can be made to the proposed PPPM. higherorder finite element solvers for the PM part and dipoles instead of monopole forghost boundary would improve the accuracy further, potentially with a smallergrid making for even faster performance. For very large problems, with billions ofparticles, there may be interesting wrinkles in making a distributed version.The PPPM philosophy could also be extended to higher-order boundary in-tegrals, or diffusion kernels, potentially to accelerate applications in and outsidefluids: preconditioning or posing sub-domain B.C. in large scale domain decom-position solvers, extending Eulerian simulation to open space, fast image and ge-ometry processing techniques.72Chapter 5ConclusionVelocity-pressure and vorticity-velocity equations are different yet equivalent formsof the Navier-Stokes equations, and it is true that with adequately accurate numer-ical methods, no significant difference should be observed when solving either ofthe two forms. However, in practice we have observed that the vorticity-velocityapproach is much more successful in cheaply resolving vorticity-rich phenomena,while the velocity-pressure formulation lends itself more easily to resolving bound-ary layers. Classifying and tracking physical phenomena with different characteris-tics using different formulations promises to be an efficient way to obtain solutionsnot possible by dealing with only one form.5.1 SummaryIn Chapter 2, a pure Eulerian solver was used to both resolve the boundary layervelocity field and the wake. In such flows what happens at the boundary domi-nates the dynamics off the boundary. The proposed method captures such flowsaccurately and cheaply.In Chapter 3, we augment velocity-pressure solvers with a correction from thevorticity-velocity formulation. This easy-to-integrate procedure is able to enhancea wide variety of existing Eulerian solvers in computer graphics, enabling themto capture vortex structures with coarser grids and much larger time steps thanpossible with other known solvers.73In Chapter 4, a pure vorticity-velocity solver is demonstrated with a PPPMscheme to reduce the N-Body summation from O(N2) to O(N). Besides the com-putational power of the PPPM method, its ease of implementation also stands outamong existing sophisticated fast summation schemes.5.2 Future workThe question remains if there are more convenient and efficient ways to bring theadvantages of vorticity-velocity methods to the velocity-pressure schemes favouredin computer graphics. IVOCK demonstrates there is potential for such schemes,but was prone to stability problems near solid boundaries at high Reynolds numberwhich weren’t satisfactorily resolved.One of the advantages of Lagrangian vortex methods is efficiently capturingflows in which vorticity is sparse, i.e. zero almost everywhere except on a lowerdimensional set. This suggests that even for Eulerian velocity-pressure schemes,adaptive grids where vorticity drives grid refinement may be attractive.A vorticity-velocity method could also be combined with a velocity-pressurescheme to resolve complex flows in a large scale, each handling the part of flow towhich they are best suited. We have already achieved some results in this direction,as shown in Fig. 5.1. In this simulation, the vFLIP solver was used to handle com-plex flows near the object boundaries, which is cheaply and weakly coupled witha Vortex-In-Cell (VIC) solver for the off-boundary flows in open space. The VICcomputation was constructed on sparse grid cells reflecting the sparse nature of vor-ticity, constructed as a collection of small fixed-sized cubic tiles of voxels, takenfrom a conceptually infinite grid but only instantiated where vorticity is present.This hybrid solver is by far the most efficient solver considered in the course ofthis thesis for capturing small-scale dynamics in a large domain, with minimumloss of detail.Finally, while free surface boundary conditions are especially difficult to in-corporate into vortex methods, compared to velocity-pressure schemes, there aremany cases in water simulation where vorticity is concentrated only on the sur-face. Potential flow schemes in particular have been extremely successful for mod-eling moderate ocean waves, while falling short of being able to properly handle74Figure 5.1: 3D simulation of flow past sphere at Re=8000. Top: one framefrom the simulated result. Middle: zoom-in of the wake pattern. Bot-tom: near-boundary flow is captured accurately with the vFLIP solverand seamlessly coupled with the free-space solution.75overturning waves and more involved solid-interaction with splashing etc. We sus-pect hybrid solvers where a sparse vortex approach handles the deep water and avelocity-pressure method handles the free surface, may provide a solution.76Bibliography[1] C. R. Anderson. A method of local corrections for computing the velocityfield due to a distribution of vortex blobs. J. Comp. Physics, 62:111–123,1986. → pages 47[2] R. Ando, N. Thürey, and C. Wojtan. Highly adaptive liquid simulations ontetrahedral meshes. ACM Trans. Graph. (Proc. SIGGRAPH 2013), July2013. → pages 5[3] R. Ando, N. Thürey, and C. Wojtan. A dimension-reduced pressure solverfor liquid simulations. EUROGRAPHICS 2015, 2015. → pages 9[4] A. Angelidis and F. Neyret. Simulation of smoke based on vortex filamentprimitives. In Symposium on Computer Animation, pages 87–96, 2005. →pages 48[5] J. Barnes and P. Hut. A hierarchical O(N log N) force-calculation algorithm.Nature, 324:446–449, 1986. → pages x, 32[6] C. Batty, F. Bertails, and R. Bridson. A fast variational framework foraccurate solid-fluid coupling. ACM Trans. Graph. (Proc. SIGGRAPH), 26(3):100, 2007. → pages 5, 8, 13, 39[7] R. Bridson. Fluid Simulation for Computer Graphics. A K Peters / CRCPress, 2008. → pages 4, 24[8] T. Brochu, T. Keeler, and R. Bridson. Linear-time smoke animation withvortex sheet meshes. In Proc. ACM SIGGRAPH / Eurographics Symp.Comp. Animation, pages 87–95, 2012. → pages 26, 27, 49, 64[9] A. J. Chorin. Numerical study of slightly viscous flow. Journal of FluidMechanics Digital Archive, 57(04):785–796, 1973. → pages 5, 45, 64, 6677[10] J. M. Cohen, S. Tariq, and S. Green. Interactive fluid-particle simulationusing translating Eulerian grids. In ACM Symp. I3D, pages 15–22, 2010. →pages 53[11] G.-H. Cottet and P. Koumoutsakos. Vortex methods - theory and practice.Cambridge University Press, 2000. → pages 48, 60[12] B. Couet, O. Buneman, and A. Leonard. Simulation of three-dimensionalincompressible flows with a vortex-in-cell method. Journal ofComputational Physics, 39(2):305 – 328, 1981. ISSN 0021-9991. → pages2, 46[13] E. Edwards and R. Bridson. Detailed water with coarse grids: combiningsurface meshes and adaptive discontinuous Galerkin. ACM Trans. Graph.(Proc. SIGGRAPH), 33(4):136:1–9, 2014. → pages 8, 25[14] S. Elcott, Y. Tong, E. Kanso, P. Schröder, and M. Desbrun. Stable,circulation-preserving, simplicial fluids. ACM Trans. Graph., 26(1), Jan.2007. ISSN 0730-0301. → pages 42[15] R. E. English, L. Qiu, Y. Yu, and R. Fedkiw. Chimera grids for watersimulation. In Proceedings of the 12th ACM SIGGRAPH/EurographicsSymposium on Computer Animation, SCA ’13, pages 85–94, New York, NY,USA, 2013. ACM. ISBN 978-1-4503-2132-7. → pages 8[16] R. Fedkiw, J. Stam, and H. W. Jensen. Visual simulation of smoke. In Proc.ACM SIGGRAPH, pages 15–22, 2001. → pages 9, 24, 27, 46[17] B. E. Feldman, J. F. O’Brien, and O. Arikan. Animating suspended particleexplosions. ACM Trans. Graph. (Proc. SIGGRAPH), 22(3):708–715, 2003.→ pages 27, 42[18] F. Ferstl, R. Westermann, and C. Dick. Large-scale liquid simulation onadaptive hexahedral grids. Visualization and Computer Graphics, IEEETransactions on, PP(99):1–1, 2014. ISSN 1077-2626. → pages 9, 49[19] N. Foster and R. Fedkiw. Practical animation of liquids. In Proc. ACMSIGGRAPH, pages 23–30, 2001. → pages 24, 27[20] M. N. Gamito, P. F. Lopes, and M. R. Gomes. Two-dimensional simulationof gaseous phenomena using vortex particles. In Computer Animation andSimulation ’95, Eurographics, pages 3–15. Springer Vienna, 1995. → pages578[21] A. Golas, R. Narain, J. Sewall, P. Krajcevski, P. Dubey, and M. Lin.Large-scale fluid simulation using velocity-vorticity domain decomposition.ACM Transactions on Graphics (TOG), 31(6):148, 2012. → pages 5, 8[22] S. Green. Cuda particles. nVidia Whitepaper, 2(3.2):1, 2008. → pages 53[23] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J.Comp. Physics, 73:325–348, 1987. → pages 46[24] N. A. Gumerov and R. Duraiswami. Fast multipole methods on graphicsprocessors. Journal of Computational Physics, 227:8290 – 8313,2008/09/10/ 2008. → pages 49[25] F. H. Harlow and J. E. Welch. Numerical calculation of time-dependentviscous incompressible flow of fluid with free surfaces. Physics of fluids, 8:2182–2189, 1965. → pages 31[26] R. D. Henderson. Scalable fluid simulation in linear time on shared memorymultiprocessors. In Proceedings of the Digital Production Symposium,DigiPro ’12, pages 43–52, New York, NY, USA, 2012. ACM. ISBN978-1-4503-1649-1. → pages 49[27] R. W. Hockney and J. W. Eastwood. Computer Simulation Using Particles.Taylor & Francis, Jan. 1989. ISBN 0852743920. → pages 45[28] A. Jacobson, L. Kavan, , and O. Sorkine-Hornung. Robust inside-outsidesegmentation using generalized winding numbers. ACM Trans. Graph., 32(4):33:1–33:12, 2013. → pages 46[29] C. Jiang, C. Schroeder, A. Selle, J. Teran, and A. Stomakhin. The affineparticle-in-cell method. ACM Trans. Graph., 34(4):51:1–51:10, July 2015.ISSN 0730-0301. → pages 8, 22[30] M. M. Kazhdan, M. Bolitho, and H. Hoppe. Poisson surface reconstruction.In Symp. Geometry Processing, pages 61–70, 2006. → pages 46, 49, 69[31] B. Kim, Y. Liu, I. Llamas, and J. Rossignac. FlowFixer: Using BFECC forfluid simulation. In Proc. First Eurographics Conf. on Natural Phenomena,NPH’05, pages 51–56, 2005. → pages 1, 23, 27[32] T. Kim, N. Thurey, D. James, and M. H. Gross. Wavelet turbulence for fluidsimulation. ACM Trans. Graph. (Proc. SIGGRAPH), 27(3):50, 2008. →pages 9, 29, 4379[33] M. Lentine, W. Zheng, and R. Fedkiw. A novel algorithm for incompressibleflow using only a coarse grid projection. ACM Trans. Graph., 29(4):114:1–114:9, July 2010. ISSN 0730-0301. → pages 9[34] M. Lentine, M. Aanjaneya, and R. Fedkiw. Mass and momentumconservation for fluid simulation. In Proc. ACM SIGGRAPH / EurographicsSymp. Comp. Anim., pages 91–100, 2011. → pages 1, 25[35] J. Lienhard. Synopsis of Lift, Drag, and Vortex Frequency Data for RigidCircular Cylinders. Bulletin (Washington State University. College ofEngineering. Research Division). Technical Extension Service, WashingtonState University, 1966. → pages ix, 17[36] C. H. Liu and D. J. Doorly. Vortex particle-in-cell method forthree-dimensional viscous unbounded flow computations. InternationalJournal for Numerical Methods in Fluids, 32(1):23–42, 2000. ISSN1097-0363. → pages 34[37] F. Losasso, F. Gibou, and R. Fedkiw. Simulating water and smoke with anoctree data structure. ACM Trans. Graph. (Proc. SIGGRAPH), 23(3):457–462, 2004. → pages 5, 27, 71[38] M. Macklin and M. Müller. Position based fluids. ACM Trans. Graph.(Proc. SIGGRAPH), 32(4):104, 2013. → pages 42[39] S. Mas-Gallic. Particle approximation of a linear convection-diffusionproblem with neumann boundary conditions. SIAM J. Numer. Anal., 32(4):1098–1125, Aug. 1995. ISSN 0036-1429. → pages 6[40] A. McAdams, E. Sifakis, and J. Teran. A parallel multigrid poisson solverfor fluids simulation on large grids. In Proc. ACM SIGGRAPH /Eurographics Symp. Comp. Anim., pages 65–74, 2010. → pages 9, 14, 15,27, 49[41] J. Molemaker, J. M. Cohen, S. Patel, and J. Noh. Low viscosity flowsimulations for animation. In SCA ’08: Proceedings of the 2007 ACMSIGGRAPH/Eurographics symposium on Computer animation, pages 9–18.Eurographics Association, 2008. → pages 8, 49[42] P. Mullen, K. Crane, D. Pavlov, Y. Tong, and M. Desbrun.Energy-preserving integrators for fluid animation. ACM Trans. Graph.(Proc. SIGGRAPH), 28(3):38:1–38:8, 2009. → pages 25, 4280[43] D. Q. Nguyen, R. Fedkiw, and H. W. Jensen. Physically based modeling andanimation of fire. ACM Trans. Graph. (Proc. SIGGRAPH), 21(3):721–728,2002. → pages 24, 27[44] S. I. Park and M.-J. Kim. Vortex fluid for gaseous phenomena. In Proc.ACM SIGGRAPH / Eurographics Symp. Comp. Animation, pages 261–270,2005. → pages 2, 7, 26, 27, 46[45] T. Pfaff, N. Thuerey, A. Selle, and M. Gross. Synthetic turbulence usingartificial boundary layers. ACM Trans. Graph., 28(5):121, 2009. → pages 7,19, 43, 48[46] T. Pfaff, N. Thuerey, and M. Gross. Lagrangian vortex sheets for animatingfluids. ACM Trans. Graph. (Proc. SIGGRAPH), 31(4):112:1–8, 2012. →pages 27, 49[47] A. Ralston. Runge-Kutta methods with minimum error bounds.Mathematics of Computation, 16(80):431–437, 1962. → pages 12[48] H. Schechter and R. Bridson. Evolving sub-grid turbulence for smokeanimation. In Proceedings of the 2008 ACM/Eurographics Symposium onComputer Animation, 2008. → pages 9[49] A. Selle, N. Rasmussen, and R. Fedkiw. A vortex particle method forsmoke, water and explosions. ACM Trans. Graph. (Proc. SIGGRAPH), 24(3):910–914, 2005. → pages 28, 48[50] A. Selle, R. Fedkiw, B. Kim, Y. Liu, and J. Rossignac. An UnconditionallyStable MacCormack Method. J. Scientific Computing, 35(2–3):350–371,2008. → pages 1, 23, 25, 27[51] R. Setaluri, M. Aanjaneya, S. Bauer, and E. Sifakis. SPGrid: A sparse pagedgrid structure applied to adaptive smoke simulation. ACM Trans. Graph.(Proc. SIGGRAPH Asia), 33(6):205:1–205:12, 2014. → pages 5, 14, 27[52] O.-y. Song, D. Kim, and H.-S. Ko. Derivative particles for simulatingdetailed movements of fluids. IEEE Trans. Vis. Comp. Graph., 13(4):711–719, 2007. → pages 8[53] J. Stam. Stable fluids. In Proc. ACM SIGGRAPH, pages 121–128, 1999. →pages 1, 23, 24, 26, 46[54] J. Steinhoff and D. Underhill. Modification of the Euler equations for“vorticity confinement”: Application to the computation of interactingvortex rings. Physics of Fluids, 6:2738–2744, 1994. → pages 9, 2881[55] M. J. Stock and A. Gharakhani. Toward efficient GPU-accelerated N-bodysimulations. AIAA paper, 608:7–10, 2008. → pages 2[56] M. J. Stock and A. Gharakhani. A GPU-accelerated boundary elementmethod and vortex particle method. In AIAA 40th Fluid DynamicsConference and Exhibit (July 2010), page 1, 2010. → pages 5, 7[57] M. J. Stock, A. Gharakhani, and C. P. Stone. Modeling rotor wakes with ahybrid OVERFLOW-vortex method on a GPU cluster. In 28th AIAA AppliedAerodynamics Conference, volume 20, 2010. → pages 2, 8[58] E. G. Tabak. Vortex stretching in incompressible and compressible fluids.Courant Institute, Lecture Notes (Fluid Dynamics II), 2002. → pages 40[59] M. Teschner, B. Heidelberger, M. Mueller, D. Pomeranets, and M. Gross.Optimized spatial hashing for collision detection of deformable objects. InProc. VMV, pages 47–54, 2003. → pages 53[60] T. Theuns. Parallel PPPM with exact calculation of short range forces.Computer Physics Commun., 78(3):238 – 246, 1994. ISSN 0010-4655. →pages 53, 59, 60[61] M. Van Dyke. An album of fluid motion. Parabolic Press, 1982. → pages 17[62] J. H. Walther. An influence matrix particle-particle particle-mesh algorithmwith exact particle-particle correction. J. Comp. Physics, 184:670–678,2003. → pages 54[63] H. Wang, K. A. Sidorov, P. Sandilands, and T. Komura. Harmonicparameterization by electrostatics. ACM Trans. Graph., 32(5):155:1–155:12,Oct. 2013. ISSN 0730-0301. → pages 45[64] S. Weißmann and U. Pinkall. Filament-based smoke with vortex sheddingand variational reconnection. ACM Trans. Graph. (Proc. SIGGRAPH), 29(4):115, 2010. → pages 26, 27, 48[65] L. Yaeger, C. Upson, and R. Myers. Combining physical and visualsimulation–creation of the planet Jupiter for the film “2010”. Proc.SIGGRAPH, 20(4):85–93, Aug. 1986. ISSN 0097-8930. → pages 2, 26, 48[66] X. Zhang and R. Bridson. A PPPM fast summation method for fluids andbeyond. ACM Trans. Graph. (Proc. SIGGRAPH Asia), 33(6):206:1–11,2014. → pages 5, 2782[67] X. Zhang, R. Bridson, and C. Greif. Restoring the missing vorticity inadvection-projection fluid solvers. ACM Trans. Graph., 34(4):52:1–52:8,July 2015. ISSN 0730-0301. → pages 22[68] Y. Zhu and R. Bridson. Animating sand as a fluid. ACM Trans. Graph.(Proc. SIGGRAPH), 24(3):965–972, 2005. → pages 1, 5, 8, 23, 25, 27, 6483"@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2017-09"@en ; edm:isShownAt "10.14288/1.0347254"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Computer Science"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivatives 4.0 International"@* ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/4.0/"@* ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Tackling small-scale fluid features in large domain for computer graphics"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/61435"@en .