Simulating SPH Fluid with Multi-Level Vorticity

The vortical motion of fluid in various scales is one of the most important elements in capturing the vivid realism of turbulent water. However, it is still challenging to resolve multi-level vorticity effectively. This paper presents a novel hybrid method that combines a smoothed particle hydrodynamics (SPH) system with multiple Eulerian grids to reproduce the multi-level vorticity. In our hybrid method, the SPH system is responsible for resolving flow velocity while the multiple grids support the SPH system in efficiently detecting and computing the multi-level vorticity field. Index Terms — Fluid simulation; Multi-level approach; Vorticity; SPH;


INTRODUCTION
The vortical motion of fluid is frequently observed both in real life and in computer animation. Examples include flow from a faucet, pouring water into a glass, and the natural stream of a river. With the vorticity, the flows often become turbulent and expose chaotic vortical motions. To capture this visually interesting effect, it is very important to resolve the vorticity in various scales. However, as the turbulence has a wide spectrum of velocity scales, it is generally difficult to capture the vivid reality through a numerical simulation. For example, high resolution is required for the spatial discretization which incurs a high computational cost.
The cost can be effectively reduced by using an adaptive approach to the simulation. However, a system using an adaptive method is still insufficient to simulate the multi-level vorticity [1] when it relies solely on the Eulerian approach or the Lagrangian approach. For example, Eulerian grid-based systems often employ more accurate numerical schemes or exploit an adaptive grid structure to overcome the difficulty in simulating the vortical motion. However, the required scales to resolve the small vortical motion are often smaller than the cell size of an adaptive grid, and the effects of the numerical diffusion are not negligible for the small scales. Similarly, the Lagrangian particle system selectively increases the number of sampling particles at the regions of in terest, such as near the free surface or around objects. However, the previous methods [2,3,4,5] emphasized only the spatial adaptiveness and paid little attention to the effects that last for a relatively long period of time, such as large scale eddies typically observed in water streams. It is inherently difficult to define time clustering operators for spatio-temporal adaptiveness due to the Lagrangian property of particles. The limitation makes the previous adaptive methods unsuitable for simulating vortical motions of fluid that vary both in spatial and temporal scales. Moreover, in spite of using an adaptive sampling method, the effects of artificial viscosity that occur during interpolation operations still exist in a smoothed particle hydrodynamics (SPH) based system. In addition, Lagrangian particle methods suffer from a difficulty in determining accurately where and how much the vorticity should be enhanced or preserved [6]. For instance, in the typical vortex particle method [7], particles with a random initial vorticity are seeded at user specified regions in a certain frame. Without a spatial constraint, the vortex particles can freely move around in the whole domain, which often results in deformation at the unwanted location. As movement of the vortex particles are determined by Lagrangian advection and cannot be predicted before running the simulation, the deformation results are sometimes unpredictable and can look unnatural in some regions.
As an alternative, a hybrid approach can be employed to resolve the multi-level vorticity. The advantages from the Eulerian and Lagrangian methods can be combined to complement each other. For example, the vorticity in a scale that is smaller than the Eulerian grid size can be resolved by Lagrangian particles, which can diminish the effect of numerical diffusion from the grid system. The location and magnitude of the vorticity can easily be determined by taking the curl to the velocity field on the grid.
Based on the hybrid idea, this paper introduces a novel method that efficiently captures the multi-scale vortical motion in SPH fluid. A SPH system combined with multi-level Eulerian grids is used to compute the vorticity in various scales. Compared to a typical hybrid system that exploits an Eulerian grid system to handle large bodies of water but uses a Lagrangian particle system to generate small scale details, such as underwater bubbles [8] or water droplets in the air, our approach utilizes the particle system to resolve fluid velocities in every scale and exploits multiple Eulerian grids to support the particle The system overview. a) The particle velocities from a SPH system are transferred to the cells in multiple grids. b) Using the cell velocity, vorticity is computed at each cell. c) The vorticity fields from each level are combined and transferred to the particles in the SPH system. d) The multi-level particle vorticity is enhanced by applying the vorticity confinement to each particle.
system to efficiently resolve multi-level vorticity. In our setting, the multi-level vorticity is computed at each level by exploiting the regularity of an Eulerian grid as illustrated in Fig. 1. The detailed vortical motion is then reproduced by enhancing the vorticity around particles using a particle-based vorticity confinement method similar to [8].
We show examples with vortical motion in highly deformable regions to demonstrate that our system effectively reproduces multi-level vorticity for SPH fluid.

SPH
The original SPH [9] is an interpolation method that solves an astrophysics problem. The method was introduced to computer graphics by Desbrun and Gascuel [10] and employed in various applications including water simulation [11], multiple fluids [12,13], viscoelastic fluid [14] and frothing liquid [15]. Despite its versatility, SPH suffers from the difficulty of enforcing incompressibility due to the inherent property of compressibility. Several methods [16,17,18] have been proposed to solve the problem. Recently, parallel and real-time simulations of SPH using GPU [4,5,19,20] have gained popularity by exploiting the locality and the simplicity of the algorithm. Our particle system is also implemented with CUDA [21].

Vortex particle
The vortex particle method was introduced to the graphics community by Gamito et al. [22] and has received much attention for simulating turbulent motion. For example, a grid-particle hybrid method for water and smoke turbulence was proposed in [7]. Similarly, a Lagrangian-based method [23] was proposed to solve the vortex formulation including boundary conditions. To reproduce realistic turbulence effects around obstacles, vortex particles with energy dynamics equations were exploited at the precomputed boundary regions [24]. Recently, hybrid methods [25,26] were proposed that exploit the vortex particles to synthesize small scale turbulence details on a low-resolution grid.

Adaptive approach
The basic idea of the adaptive approach is to change the particle resolution for the region of interest by spatial sampling. In computer graphics, Adams et al. [2] used the approach in water simulation by changing the particle resolution according to the geometric feature size. For example, a high resolution creates small splashes or thin sheet effects while a low resolution creates internal water volumes. Techniques with elaborate sampling methods, such as using a high pressure gradient [3] or a local pressure ratio and the number of neighbors [4,5] were recently proposed.

Hybrid approach
A hybrid approach combines an Eulerian grid system and a Lagrangian particle system together to complement each other. Several variants have been proposed differing the degree of importance of each component to achieve a specific purpose. Early methods utilize auxiliary particles to assist the grid system with surface tracking or with correction for interface capturing [27]. In an effort to reproduce small scale details, use of particles started to replace the dissipative grid in the advection process [28]. Escaped particles are incorporated into the SPH system to simulate underwater bubbles [8]. A two-way coupling method [29] was proposed to simulate both dense and diffuse fluid regions using the grid and SPH system, respectively.
Our system is a variant of the particle-grid approach, which is the opposite of previous grid-particle approaches in a hybrid domain. The particle-grid approach is mainly based on a particle system, which has the advantage of preserving sub-grid scale features from numerical diffusion. In contrast, the grid-particle methods require an expensive grid refinement process or an additional reseeding process to capture the sub-grid scale details. Another advantage of the particle-grid approach is that most of the simulation process can be parallelized easily, and can gain speedup through hardware acceleration. We focus on building a simple but effective solution that takes full advantage of the modern particle system while maintaining simplicity in the hybrid process. The steps for coupling are simplified and the parallel processing is maximally exploited.

Incompressible SPH
The behavior of the SPH fluid is described by the approximation of the momentum equation where i x is the particle position, i v is the particle velocity, p is the pressure, µ is the viscosity constant, and ext i f is an external force term for gravity, user defined forces, or vorticity confinement forces.
i ρ , p ∇ and v ∆ symbolize the interpolation kernel-based approximation of the density field, pressure force field, and viscous force field, respectively. In our implementation, we employ a poly kernel for field interpolation, a spiky kernel to approximate the gradient, and a viscosity kernel for the computation of the Laplacian, respectively. In order to minimize the effects of spurious bouncing motion during the particle simulation, we apply the predictive-corrective incompressible SPH (PCISPH) solver [18]. In PCISPH, the pressure value for each particle is updated by solving the following prediction-correction scheme iteratively.
ρ is the predicted density, δ is the scaling variable, and 0 ρ is the rest density. In our test scenes, we precompute the δ value using a randomly picked internal particle. With the fixed scaling factor δ , we repeat the prediction-correction step at each frame with 3 iterations to enforce the incompressibility of the SPH fluid as in [18].

Multi-level grids
To resolve the multi-level vortical motion stably and efficiently, we exploit the spatially uniform Eulerian MAC grid. We generate multiple grids with different resolutions. The number of levels n and the ratio r between each level are provided by the user. The generated grids correspond to multiple spatial sub-samplings of the domain. In the highest resolution grid, the size of a cell is the same as the SPH particle's smoothing radius and the time interval is the same as the time-step of the SPH simulation. A cell size i d and a time interval i t at each grid are determined by applying the Courant-Friedrichs-Lewy (CFL) condition to each level where i j n < ≤ and u represents the velocity, and cfl k is a system parameter. The multiple cell sizes and time intervals are required to capture the vortical motion with different scales and speeds. For example, the motion of a large eddy that occupies many cells can be captured in the low resolution grid, which has larger cells and a longer time interval.

Hybrid deformation detection
The goal of our hybrid system is to exploit the regularity of the Eulerian grid to complement the SPH particle system in resolving the multi-level vorticity. For the detection of the deformation, we propose a hybrid technique. We formulate the particle deformation detection process as a local eigen-analysis on grid cells. The Principle Component Analysis (PCA) applied to the particles in each cell computes the variation of particles with respect to the X,Y, and Z axis of the cell. Assuming that a cell with coordinate ( , , ) i j k has the center position l c in level l and number n of particles p , the covariance matrix l Cov for the cell is It is possible to encode the deformation for each particle. For example, we can change the l c to a particle centric criterion like a centroid ∑ . However, the particle-based deformation encoding often requires attention in handling an exceptional case, such as a lump of particles or outliers. While the cell based particle deformation encoding is an approximation based on lower resolution, it handles the irregular situations including deficiency of neighboring particles better than the particle-based deformation encoding. In addition, other particle centric criteria, such as the local pressure ratio, the number of neighboring particles [4], or the combination of both [5] are dependent on the accuracy of kernel interpolation, which is vulnerable to the deformation. Moreover, for different scene settings, a time consuming trial and error process is required to determine a reasonable threshold for a single criterion or optimal weights for the criteria combination.
To detect the deformable regions, we examine the eigenvalues of the For example, if the particles are in rapid divergent or convergent motion, the change of variance in each direction will be high enough to mark the container cell as deformable.

Hybrid vorticity enhancement
Once deformable cells are detected, we can compute the multi-level vorticity using multiple grids. The process begins by constructing a velocity field for each grid. A weighted average of the particle velocities produces the cell velocity u .
where i v is the particle velocity and i d is the distance be-

∑ ∑
By accumulating the multiple vorticity fields, we get a single vorticity field $\omega$. The particle vorticity l ω is computed by tri-linear interpolation on the field. For each particle, the vorticity confinement force is computed only using neighboring particles. Our system was implemented using CUDA to maximize the usage of parallel computation. The most frequent and time-consuming process in the particle system is to search for the nearest neighbors. A particle hashing is applied to the auxiliary grids followed by the radix sort algorithm for efficiency. We can achieve a real-time performance on the simulation for our test examples. We also use a CUDA-based marching cube with a grid (128x128x128) to gain a speed-up in the mesh generation process, which now takes only several seconds in our system.
We reuse the grids to calculate the cell-based velocity and the vorticity. The multi-level grids are extended auxiliary grids with different resolutions. They are stored in the GPU memory as 1D array. An individual CUDA thread is activated and run in parallel for each cell's coupling process. For example, the vorticity confinement in SPH evokes as many CUDA threads as the number of cells to utilize the velocity and the vorticity defined in each cell. The association of the cell and the particles inside is registered in the GPU memory and updated at every frame. The usage of the GPU accelerates the particle-grid coupling process. In addition, we further improve the performance when solving for the eigenvalue of a 3-by-3 matrix by implementing a simple kernel function in the GPU.
The second column in Fig. 4 shows the result of a simulation of a single dam-breaking scene. The simulation was carried out using 44,649 particles. We can notice that the surface has more deformation with our method. Although, the vortical motions are somewhat exaggerated by adjusting the user parameter, the effects of the multi-level vorticity make the flow look like a natural turbulence. A user can adjust the epsilon values to emphasize the vorticity at a specific level or to modify the effective ranges of the vorticity. We use 1.0, 0.15 and 0.01 as the epsilon values for the scene. While a tedious trial and error procedure is often required to determine the epsilon values, they can be reused for several scenes once they are set.
In the third column, we demonstrate a scene of pouring water on a bunny object. In contrast to the previous example, we can see the effects of the multi-level vorticity from the particles at the thin surface. After hitting the bunny, the motions of the particles become a little bit more scattered with the generated multi-level vorticity.

V. CONCLUSION
We present a new particle-grid hybrid system to efficiently reproduce multi-level vortical motion in SPH fluid. By combining a SPH system with multiple Eulerian grids, we can detect the deformable region robustly while efficiently computing the vorticity for particles. With the flexibility and the simplicity of our multiple grid system, one can easily extend our system to a wider spectrum in aspects of space in time. In a future study, we plan to extend our framework to resolve the sub-particle scale vortical motions as well..

ACKNOWLEDGEMENT
This work was supported by culture contents industry research and development program of KOCCA/MCST (210-7602-003-10743-01-007, Software Development for 2D to 3D Stereoscopic Image Conversion). tute service) and participated in several commercially shipped games including Pang-Ya online, Trickster online, and whiteday as a game programmer. As a graduate student, he participated in several research projects that are funded by government while experiencing various production -level developments including plug-in programming for Maya, Motion Builder, Real-flow and Houdini using C++. He also has more than 1 year of experience on CUDA programming. Current research interest is an efficient simulation of multi-level vorticity. As in his previous research, "Multi-level vorticity confinement for water turbulence simulation", he pursues to find practical technique to be used in real-world production pipeline with cost efficiency.