Abstract
We present a novel Smoothed Particle Hydrodynamics (SPH) based algorithm for efficiently simulating compressible and weakly compressible particle fluids. Prior particle-based methods simulate all fluid particles; however, in many cases some particles appearing to be at rest can be safely ignored without notably affecting the fluid flow behavior. To identify these particles, a novel sleepy strategy is introduced. By utilizing this strategy, only a portion of the fluid particles requires computational resources; thus an obvious performance gain can be achieved. In addition, in order to resolve unphysical clumping issue due to tensile instability in SPH based methods, a new artificial repulsive force is provided. We demonstrate that our approach can be easily integrated with existing SPH based methods to improve the efficiency without sacrificing visual quality.
1. Introduction
With an observable development of physically based methods, fluid simulations have popularized themselves in many graphics applications such as feature films and computer games. There are various physics-based fluid simulation techniques available in the context of computer graphics, but here we focus on SPH based methods due to their simplicity and flexibility.
In SPH based methods, the simulation and visual quality are largely influenced by the resolution of the fluid, that is, the number of particles. Nevertheless, the higher the resolution is, the more computing power is required. In general, only a few milliseconds can be used for physics computation in computer games. As a consequence, low resolution fluid simulation is used in games to get acceptable frame rates, for example, in NVIDIA PhysX engine [1] SPH fluid with less than 60 k particles is typically used for real time applications. In SPH based methods, regardless of its usage in off-line simulations or real time 3D games, reducing the total number of fluid particles can lead to a lower computational overhead. To this end, adaptively sampled SPH method [2] has been introduced to allocate more computing resources to specific regions such as the visually important fluid surface. Similarly, our goal is to develop an efficient fluid simulation method based on Lagrangian SPH model, which is able to lower physical simulation times without compromising the resulting visual quality.
Inspired by the sleepy strategy frequently used in rigid body simulations, we apply this idea in the context of SPH based methods. However, unlike adaptive sampling techniques which also concentrate computational resources only on visually interesting regions, we only use a single resolution. We begin by identifying different particle types according to physical or geometrical criteria. Next, we perform only parts of the whole physical simulation pipeline for different particle types. Although we show that our model can reduce computational overhead compared to standard WCSPH approach while still preserving plausible visual results, our algorithm can be also integrated with standard SPH method commonly applied in real time physics engines.
Additionally, we note that traditional SPH based methods usually lead to particle clumping as a result of density underestimation near the fluid surface where negative pressures are generated that causes the particles to cluster, which is known to be tensile instability issue. We avoid this problem by introducing a positive repulsive force to naturally separate them apart. The experiment results show that our repulsive force can eliminate particle clumping issue that occurred on the fluid surface in previous methods.
The main contributions of this work are as follows:(1)an efficient SPH based algorithm grounded on a novel sleepy strategy, which can significantly enhance the performance of fluid simulation;(2)a new artificial repulsive force to resolve tensile instability issues on account of inaccurate density estimation, which can nicely avoid particle clumping near the fluid surface.
In the remainder of this paper, we first briefly review the related work in Section 2. In Section 3, we give an overview of the standard SPH and WCSPH framework. Section 4 continues with a discussion of our sleepy method and a new artificial repulsive force for resolving tensile instability issue in SPH fluids. Implementation details, experimental results, and analysis are given in Section 5. Finally, conclusion of the work is presented in Section 6.
2. Related Work
SPH was first independently developed by Monaghan and Gingold [3] and Lucy [4] to solve astrophysical problems. It was introduced to the graphics community by Stam and Fiume [5] to simulate fire. Muller et al. [6] used this method to interactively simulate compressible fluids. Since then, various fluid effects have been simulated using this technique in computer graphics, such as fluid-deformable interaction [7, 8], fluid-fluid interaction [9, 10], fluid-solid coupling [11, 12], melting [13, 14], surface tension and cohesion [15, 16], weakly compressible fluids [17, 18], and incompressible fluids [19–22].
Most of the work mentioned above adopts a single resolution model both in terms of space and time, allocating computational resources uniformly to each particle which is not always necessary. In order to simulate large-scale fluid scenarios in reasonable time on a desktop computer, several adaptive sampling techniques have been introduced.
One way that employs the concept of adaptive sampling is adaptive particle refinement, which devotes most of the computational resources to visually or physically important fluid areas. Desbrun and Cani [23] used an adaptive scheme in the context of SPH where particle sampling and time sampling are adapted according to a given compromise between accuracy and efficiency. Adams et al. [2] have presented a spatial adaptive particle fluids approach in which particles are subdivided in geometrically complex and visually important areas and merged in regions of low geometric complexity and of low visual interest, and a smooth particle level distribution is also maintained to keep the physical simulation stable. This method can significantly reduce computational cost, making it possible to simulate large-scale and more complex scenarios in acceptable time. However, one drawback of this technique is that the neighbor search is a major performance bottleneck since hierarchical data structures are expensive to construct and maintain. In addition, particle splitting and merging may lead to inaccurate density and force profiles which are important for the accuracy and stability. Zhang et al. [24] extended this method to simulate weakly compressible fluids; they achieve interactive frame rates by executing both the simulation and rendering fully on the GPU. Yan et al. [25] added a physical importance criterion for adaptively splitting or merging particles and also used GPUs to accelerate the simulation. Orthmann and Kolb [26] applied temporal blending scheme to reduce the total number of particles while still being capable of using large time steps. They can simulate up to a million of particles at interactive speed by exploiting the massively parallel computation capabilities of GPUs.
Another way towards adaptive sampling is adaptive time discretization, which uses different time steps for each simulation frame. Desbrun and Cani [23] automatically adjusted the individual time step for each particle with respect to the Courant-Friedrichs-Lewy (CFL) condition; thus instability problem is avoided. Ihmsen et al. [27] have proposed an adaptive time-stepping method for predictive-corrective incompressible SPH (PCISPH), which adaptively evaluates appropriate time steps regardless of the scenario. By utilizing this temporal adaptivity, the overall computational overhead can be largely reduced compared to constant time-stepping techniques.
Recently, Solenthaler and Gross [28] used a multiscale sampling technique that uses two coexisting simulations with different resolution levels to cut down the overall computational overhead, which does not need to split or merge particles, and can stably simulate two simulations of largely different resolutions. This two-scale particle approach can improve the simulation efficiency linearly to the reduction of the total particle count while still achieving similar surface fluid flow behavior. However, the mass conservation is lost because particles need to be dynamically added or removed in the fine resolution region. Horvath and Solenthaler [29] extended this method to simulate high resolution incompressible fluids with multiple levels, which further reduces the computational costs and also conserves the total mass in high resolution area.
3. Background
SPH is a particle-based interpolation method which discretizes continuous field using the properties defined in all particles. As for SPH fluids, this interpolation can be described as The value denotes the physical property, the density, the mass, the distance between the evaluated particle and its neighbor , the weighting kernel, and the smoothing length.
In standard SPH, the pressure is computed from the density fluctuations of the fluid particles through ideal gas equation [6], which leads to perceivable compressibility artifact. To resolve this issue, WCSPH is introduced into graphics community [18]. In WCSPH, the Tait equation is applied to enforce incompressibility of the fluids. While Tait equation is essentially an equation of state (EOS), it uses large stiff values to generate huge penalty force to avoid much smaller distance between fluid particles usually appearing in the compressible case. The basic SPH/WCSPH method is summarized in Algorithm 1. At the beginning of each physics simulation step, the neighboring particles need to be found and then used to calculate particle’s physical attributes, such as density, pressure, and various forces acting on them. Particle’s density is evaluated through the following equation: The pressure is calculated from Tait equation: where is the pressure at particle , is a stiffness parameter to allow small density fluctuation, is the reference density, and is a user defined value which is usually set to 1 for SPH and 7 for WCSPH.
| ||||||||||||||||||||||||
The pressure force and viscosity force on each particle are given by the following equations: where and denote individual velocities of particle and .
4. Sleepy Method
In standard WCSPH, the physical quantities of all particles are being traced and updated at each time step to obtain a correct and stable simulation, whereas this is not always necessary in terms of visual plausibility. Considering the rigid body simulations, some rigid bodies are regarded as sleepy if they satisfy specific criteria, meaning that we do not need to compute their physical properties. Similarly, we observe that in WCSPH, there are usually many particles that can be treated as quasi static, which means their moving states can be safely ignored during several consecutive time steps; thus lots of computational resources can be saved. In the following, we will discuss the sleepy strategy used in our framework as well as the solution to tensile instability problem in the SPH based methods.
4.1. Sleepy Strategy
In adaptively sampled SPH method, particles are categorized to different levels according to their radii. Similarly, we have partitioned the whole fluid particles into three different types, that is, nonsleeping particles, intermediate particles, and sleeping particles, as shown in Figure 2. However, we still use a regular particle sampling which is better at reproducing physical properties and more efficient at neighbor searching compared to spatial adaptive sampling method. This is because that the latter uses hierarchical data structures such as kd-trees to perform neighbor search which are expensive to construct.
The processing of nonsleeping particles is the most computationally demanding part of the overall simulation. The computing demand of intermediate particles is second only to nonsleeping ones since they are needed to model the virtual boundary conditions for nonsleeping particles. The sleeping particles mark the most important ones in our method because it is precisely on the sleeping particles that a large amount of computational overhead can be saved.
Similar to its concept in rigid body simulation, a sleepy particle does not move, nor does its current state imply that it means to move in the instant future. Thus computational resources are saved by ignoring calculation of their physical properties. Intuitively, these sleepy particles should be deep inside the fluid or under a slow motion, while the nonsleeping particles normally lie in physically complex or visually important regions such as the fluid surface.
As for the heuristic to determine nonsleepy particles, we choose the following criteria: where is the kinetic energy of particle , is the local pressure, is the distance from fluid particle to the normalized center of mass of its neighbors, and , , and are user defined values. Here we adopt three criteria for detecting nonsleepy particles:(i)kinetic energy criterion means that the particles in fast moving fluid region should be marked as nonsleepy;(ii)pressure criterion means that the particles in fluid region with high turbulence should be marked as nonsleepy;(iii)distance criterion means that the particles near the fluid surface or boundary should be marked as nonsleepy.
Now that the nonsleepy particles can be identified through (6), the other two particle types can also be easily identified afterwards as will be discussed in Section 4.3
4.2. Tensile Instability
In standard SPH based methods, when fluid particles are under tensile stress state, their motions become unstable which appears in the form of either particle clumping or complete blowup. This instability is usually referred to as tensile instability. This is due to the negative pressures of the surface particles that have fewer neighbor particles than needed to maintain the rest density (see Figure 3(a)).
In SPH literature, truncating the pressure to nonnegative is a commonly used scheme to avoid tensile instability, but with the side effect that the particle cohesion will be reduced. In this paper, we introduce an artificial repulsive force to solve this problem: where is the repulsive coefficient, and are pressures of particle and , is the distance between these two particles, and is the average distance from its neighbors to particle .
When the pressures of fluid particles become negative, we use this artificial repulsive force instead of pressure force (4). Thus, unphysical attraction force as a result of the negative pressure becomes a positive repulsive force which can effectively remove the tensile instability. Figure 3 demonstrates that particle clumping problem that arises because of attraction forces between surface particles with negative pressures is avoided.
4.3. Sleepy SPH/WCSPH Algorithm
Algorithm 2 is our sleepy algorithm for compressible and weakly compressible fluid simulations. Initially, we determine a subset of the fluid as nonsleepy particles according to (6). After that, we apply a flood fill technique to these particles. Specifically, we first perform neighbor search for nonsleepy particles. Then, those neighboring particles of the nonsleepy ones with a magnitude of kinetic energy larger than a user defined threshold are flood filled to be nonsleepy. Last, we mark the remaining neighbor particles as intermediate. Thus, similar to the reactivation concept in rigid body simulation, the sleepy fluid particles could be reactivated again by simply updating their velocities. Afterwards, the sleepy particles can be easily identified by marking all the rest particles as sleepy. After all particle types have been identified, we perform neighbor search, density, and pressure computation as well as the time integration for only a part of fluid particles.
| ||||||||||||||||||||||||||||||||
We notice that in each time step, neighbor search, density, and pressure calculation are only needed for nonsleepy and intermediate particles; force calculation and time integration are only necessary for nonsleepy particles. Hence, lots of computational overhead can be saved compared to the standard SPH/WCSPH algorithm.
5. Implementation and Results
5.1. Implementation Details
In our implementation, we use cubic spline kernel as described in [30] for all calculations of physical attributes of fluid particles. The smoothing length is set to twice the rest distance between particles to make sure there are usually 30 to 40 neighbors for most particles. For the neighbor searching, a uniform grid is used to accelerate the process since it is a major performance bottleneck of the whole physics simulation. The time step is chosen to be the largest value satisfying CFL conditions. For the complex boundary handling such as letters as shown in Figure 1, we first sample the surface of static solids with boundary particles using the method from [31] and then adopt approaches from [11] to evaluate physical properties. While for simply shaped boundaries such as walls, we apply wall weight function idea from [32] since it can reduce the computational overhead compared to [11] in this case.
(a)
(b)
(a)
(b)
For our sleepy algorithm, we perform neighbor search for nonsleeping particles together with particle type determination at the beginning of each simulation step. Thus, the computation waste of another loop for type determination can be avoided. Our code is implemented in C++, and all experiments have been performed on an Intel Core i7 3.4G CPU with 8 GB of RAM. For surface reconstruction and fluid rendering, the technique from [12] is used.
5.2. Visual and Performance Results
In this section, we present several scenes simulated with our sleepy method in the context of WCSPH, which show that our algorithm can obtain similar visual results and is more efficient than standard WCSPH. Moreover, our technique can be also easily integrated with existing SPH method usually used in game physics engines.
Figure 1 shows weakly compressible water with 500 k fluid particles flow around static letters sampled with 9 k boundary particles; the corner breaking dam collides with static letters and then comes to rest. This example is created using sleepy WCSPH together with boundary handling methods from [11, 32]. Note that there is little difference between the resulting visual quality of our method and the standard WCSPH, while our approach is typically faster. In Table 1, we make a performance comparison between standard WCSPH and our sleepy method. From the table, we can note that a speedup of 3.05 can be achieved using our method compared to standard approach.
In Figure 4, we make a visual comparison between the final effects with and without our repulsive force model based on sleepy WCSPH. From the figures we can see that the tensile instability issue is avoided and surface tension like visual effects is obtained as well as the particle distribution has been improved due to the proposed repulsive force model. Therefore, the splash effects are more naturally captured with our method as can be seen from Figure 4(b).
(a) Without our repulsive force
(b) With our repulsive force
To illustrate the scaling of our method with the portion of the sleepy particles, we animate corner dam break scenes with different number of particles (see Figure 5) and compare their performance results in Table 2. As shown, when simulating scenarios with most particles under fast moving state as is the case in Figure 5(a), the speedup of our technique is small. This is because, under this circumstance, not too much computational overhead could be saved since portion of sleepy particles is relatively small. On the other hand, if the percentage of sleepy particles is larger as is the case in Figure 5(d), a higher performance boost can be obtained. On the basis of this analysis, we can draw a conclusion that our method performs better if a larger portion of the whole fluid particles are under sleeping state.
(a) 50 k
(b) 100 k
(c) 200 k
(d) 400 k
6. Conclusion and Future Work
In this paper, we have proposed an efficient method for taking advantage of the sleepy state of particles in a compressible and weakly compressible fluid simulator that allows focusing computational resources only on parts of the whole fluid volume. Experimental results show that our method can capture similar realistic fluid effects with less computational needs. Our sleepy technique is also suitable for particle-based fluids widely applied in current generation game physics engines to get a further speedup. Additionally, we have introduced a new artificial repulsive force to avoid particle clumping as a result of neighbor deficiencies near the fluid surface.
As future work, we plan to speed up our algorithm by using multicore CPUs and many-core GPUs and by integrating adaptive time step techniques.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors thank the reviewers for their helpful comments. This work is supported by National High Technology Research and Development Program of China (no. 2012AA011503) and the preresearch project (no. 51306050102).