Adv Search
Home | Accepted | Article In Press | Current Issue | Archive | Special Issues | Collections | Featured Articles | Statistics

2021,  3 (2):   129 - 141   Published Date:2021-4-20

DOI: 10.1016/j.vrih.2018.12.001
1 Introduction2 Related work2.1 SPH foundation2.2 Adaptive sampling methods2.3 Neighbor search3 Adaptive-smoothing length method3.1 Density-based adaptive smoothing length3.2 Neighbor search scheme3.3 Interacting forces4 Experimental results4.1 Number of neighbors4.2 Density field4.3 Visual effects4.4 Simulation efficiency5 Conclusion and future work


In the smoothed particle hydrodynamics (SPH) fluid simulation method, the smoothing length affects not only the process of neighbor search but also the calculation accuracy of the pressure solver. Therefore, it plays a crucial role in ensuring the accuracy and stability of SPH.
In this study, an adaptive SPH fluid simulation method with a variable smoothing length is designed. In this method, the smoothing length is adaptively adjusted according to the ratio of the particle density to the weighted average of the density of the neighboring particles. Additionally, a neighbor search scheme and kernel function scheme are designed to solve the asymmetry problems caused by the variable smoothing length.
The simulation efficiency of the proposed algorithm is comparable to that of some classical methods, and the variance of the number of neighboring particles is reduced. Thus, the visual effect is more similar to the corresponding physical reality.
The precision of the interpolation calculation performed in the SPH algorithm is improved using the adaptive-smoothing length scheme; thus, the stability of the algorithm is enhanced, and a larger timestep is possible.


1 Introduction
In the field of physics-based fluid simulation, the smoothed particle hydrodynamics (SPH)[1] method is being increasingly employed for various applications, including entertainment technology and engineering design. The SPH method, which is particle-based, separates a continuous fluid into interacting particles[2,3]. By calculating the force on each particle, the particles are continuously moved, mimicking the flowing fluid. Owing to its convenience and flexibility in complex fluid simulations, SPH has been widely used in digital entertainment and various other fields; moreover, it has already become one of the main tools for fluid simulation.
In the SPH method, the smoothing length significantly affects the simulation efficiency and accuracy. Improper setting of the smoothing length can reduce the simulation accuracy considerably. In three-dimensional (3D) simulations, employing 30-40 particles in the support domain generally yields good accuracy. In traditional SPH, the particle support domain is fixed; thus, different regions differ with regard to the number of neighboring particles in the support domain. It is difficult to obtain a consistent kernel approximation in the calculation domain, which leads to excessive particle interpolation errors, as well as the loss of minor details during the fluid simulation. Therefore, in accordance with the idea of allocating more resources to areas of visual interest, researchers investigated adaptive sampling techniques for improving the computational efficiency without compromising the visual quality[4-6].
In this paper, an adaptive smoothing length SPH (ASLSPH) fluid simulation method is proposed for reducing the calculation error caused by the fixed smoothing length of particles. By establishing the relationship between the smoothing length and the density of particles and adjusting the smoothing length adaptively according to the changes in the particle density, the method can reduce the calculation error, improve the accuracy of interpolation calculations, and achieve simulation results that are closer to the physical reality.
A limitation of our method is that the neighbor search efficiency is lower than that of the traditional method, because the smoothing lengths differ among particles. Therefore, a fast neighbor search scheme based on an auxiliary grid is proposed for increasing the simulation speed. Additionally, to solve the problem of particle interaction asymmetry caused by the variable smoothing length, a symmetric kernel function is designed, which maintains the stability of the algorithm.
The main contributions of this paper are as follows:
• A density-based adaptive smoothing length scheme to improve the calculation accuracy;
• A symmetric interpolation scheme for the kernel function to maintain the stability of the algorithm;
• A grid-based fast neighbor search scheme for particles with a variable smoothing length.
2 Related work
2.1 SPH foundation
The first particle-based method was proposed by Miller et al.[7] Among the particle-based methods, SPH[8] is the most commonly used. In 2003, Muller et al. introduced the SPH method to the field of computer graphics[1]. Since then, scientists have continuously investigated the method, which has promoted its rapid development and progress[9-12]. The incompressibility of fluids—particularly liquids—significantly affects the visual quality of the fluid animation. Thus, scientists have proposed various pressure solvers to ensure fluid incompressibility. Becker and Teschner proposed weakly compressible SPH (WCSPH) based on the Tait equation[13]. Solenthaler and Pajarola proposed predictive-corrective incompressible SPH (PCISPH) to ensure that the density is constant throughout the simulation[14]. Ihmsen et al. proposed implicit incompressible SPH (IISPH)[15]. Bender and Koschier combined two new implicit pressure solvers and proposed divergence-free SPH (DFSPH) to limit volume compression and prevent divergence in the velocity field[16,17].
2.2 Adaptive sampling methods
Several methods have been proposed by graphics researchers to achieve better simulation results and a higher efficiency. Generally, a larger number of particles corresponds to a higher calculation accuracy, but simulations with a large number of particles are time-consuming. It is difficult to simulate a large number of particles in a given time range. Therefore, the adaptive sampling technique has been investigated for improving the allocation of computational resources.
The spatial adaptive method attempts to re-sample the particles so that the non-uniform sizes of the particles can be employed to obtain animation results of higher quality. The idea of the adaptive particle size[4,5] is to dynamically adjust the radius of the particle according to the complexity and visual importance of the fluid region. To minimize the sampling error, the Lagrange multiplier[6] is typically used to minimize the local error function in computational fluid dynamics. In computer graphics, simple and more efficient sampling operations have been used, such as particle merging or splitting[18-21]. In the high-resolution region, a particle is split into a certain number of child particles, whereas in the low-resolution region, they merge into a single parent particle.
In addition to spatial adaptive sampling, timestep adaptation is helpful for improving the simulation efficiency. The adaptive-timestep method[22-24] dynamically adjusts the integral timestep in each simulation step to automatically adapt to the Courant-Friedrichs-Lewy condition, thus improving the simulation efficiency.
The aforementioned adaptive methods were mainly used to improve the calculation efficiency, whereas the adaptive smoothing length methods are used to improve the calculation accuracy and algorithm stability. Gingold and Monaghan proposed an adaptive smoothing length method for assigning independent smoothing lengths to each particle[25]. Subsequently, Nelson and Papaloizou studied the relationship between the smoothing length variation and the momentum and energy equations; further, they proposed a variable smoothing length method to reduce the energy error and improve the calculation accuracy[26]. Sigalotti et al. introduced adaptive density kernel estimation and applied an adaptive smoothing length to free surfaces[27]. Olejnik et al. used the Okubo-Weiss parameter to update the smoothing length[28]. Ren et al. proposed dual-support SPH, which allows the use of a variable smoothing length and satisfies the conservation of momentum and energy[29]. The foregoing methods improved the simulation accuracy, but the algorithms are complex, which reduces the simulation efficiency. We propose an adaptive smoothing length algorithm based on particle density variation, which can not only improve the computational accuracy but also maintain a high computational efficiency.
2.3 Neighbor search
The neighboring particles of each particle in SPH change dynamically over time. Therefore, a neighboring particle search is required in each timestep to calculate the interpolation weight of the neighboring particles. This is usually the most time-consuming step in the SPH simulation process. For an SPH simulation with a fixed smoothing length, a uniform grid[30] can be used to improve the efficiency of the neighbor search. By reordering the particles[31], the search speed can be further improved. For SPH simulations with a multiresolution particle size or a non-uniform smoothing length, hierarchical data structures such as trees[19,32,33] are typically used to conduct the neighbor search. Gong et al. proposed a parallel background grid neighbor-search method based on adaptive mesh refinement, which can significantly improve the neighbor search efficiency[34]. Additionally, because there is almost no data dependency, the SPH method can be easily mapped to a graphics processing unit (GPU)[35,36] to perform parallel computation. A GPU-based neighbor list algorithm[37-39] can be used to simulate large-scale fluid animation owing to its performance advantages.
3 Adaptive-smoothing length method
In contrast with the Eulerian grid method, in SPH, it is difficult to ensure the incompressible property of the fluid, which may lead to an uneven particle distribution in which the particles are excessively dense or sparse in some regions. The conventional SPH fluid simulation method usually adopts a fixed particle smoothing length for the neighbor search and interpolation calculation. Owing to the uneven distribution of particles, particularly at the free surface, the fixed smoothing length of particles may lead to a large difference in the number of neighbors among particles, which affects the interpolation accuracy.
In Figure 1, particles i and j are at different positions in the flow field. The particles around i are densely distributed, e.g., in the deep-water area, while the particles around j are sparsely distributed, e.g., near the fluid surface. With a fixed smoothing length h, owing to the insufficient number of neighboring particles, the interpolation accuracy for the physical quantities on particle j is reduced. To solve this problem, an ASLSPH method based on the particle density is proposed. By establishing the relationship between the smoothing length of the particles and the weighted average of the density of the neighbors of the particles, the smoothing length of the particles can be adjusted adaptively. Even in the case of a special particle distribution, e.g., a lack of neighboring particles or an uneven distribution of particles, the interpolation can maintain a high accuracy, and this yields improved interpolation accuracy and stability of the simulation.
In the proposed method, the smoothing length of each particle is a variable so that the number of neighboring particles can be kept within a stable range; thus, a high interpolation accuracy can be maintained for each particle. Additionally, to ensure that the computed forces satisfy Newton's third law, we designed a corresponding neighbor search scheme and interpolation kernel function for maintaining the stability of the algorithm. The simulation process is presented in Algorithm 1, where the improved parts (relative to the standard SPH method) are shown in bold.

Algorithm 1: ASLSPH method

1: while simulation do

2: for each particle i do

3: calculate the smoothing length of i (Section 3.1)

4: end for

5: adaptive smoothing length neighbor search (Section 3.2)

6: for each particle i do

7: calculate the density of particle i (symmetric interpolation, Section 3.3)

8: end for

9: for each particle i do

10: calculation the pressure force, viscous force and external force of particle i (symmetric interpolation, Section 3.3)

11: end for

12: for each particle i do

13: update the velocity of particle i

14: update the position of particle i

15: end for

16: end while

3.1 Density-based adaptive smoothing length
To obtain a high approximation accuracy of the interpolation kernel for each particle, the smoothing length is adjusted adaptively according to a defined scaling factor, which is controlled by the weighted average of the density of the neighboring particles. Thus, the number of neighboring particles is maintained within a stable range. Specifically, the smoothing length
h i
of particle
is defined as
h i ' = α i h ,
represents the initial radius of each particle, and
α i
is a scale factor used to adaptively control the smoothing length.
α i
is defined as follows:
α i = ρ i ¯ ρ i 1 d
represents the dimension of the simulation (2 and 3 for two-dimensional (2D) and 3D simulations, respectively),
ρ i
represents the density of particle
, and
ρ i ¯
represents the weighted average of the density of the neighboring particles.
ρ i ¯
is given as follows:
ρ i ¯ = j = 1 n ρ j i d j i j = 1 n 1 d j i
represents the number of neighboring particles for particle
ρ j i
represents the density of the
j t h
neighboring particle, and
d j i
represents the squared Euclidean distance between particle
and its
j t h
neighboring particle.
d j i
is given as follows:
d j i = x i - x j i 2 2
x j i
represents the position of the
j t h
neighboring particle.
n = 0
, particle
is considered as an isolated particle, and we set
h i = h
. If
ρ i
ρ i ¯
, the distribution of particles around particle
is sparse, and the number of neighbors is too small. If the smoothing length is fixed, an interpolation error is generated owing to the insufficient neighboring particles. However, in our method, according to Equation (2),
α i
> 1; thus, a smoothing length of >
is eventually assigned to
h i
. With a greater smoothing length, the number of neighboring particles increases accordingly, and eventually the interpolation error is reduced. In contrast, in places where the particles are densely distributed, the smoothing length
h i
decreases to a value smaller than
, thereby avoiding an excessive number of neighboring particles. As shown in Figure 1, after the smoothing length of particle
was adaptively adjusted, the number of neighboring particles was approximately identical to that for particle
. Thus, the calculation accuracy of the interpolation kernel can be improved.
3.2 Neighbor search scheme
Owing to the inconsistency in the smoothing length of the particles, a problem with our adaptive method is that the determination of the neighboring particles is asymmetric. Suppose that the smoothing length of particle
exceeds that of particle
; i.e.,
h i > h j
. When the distance
d i j
between the two particles satisfies
h j < d i j < h i
, particle
regards particle
as its neighbor, but particle
does not. This asymmetry results in a force imbalance between particles
, which is inconsistent with Newton's third law. To solve this problem, we employ the average value of
h i
h j
, i.e.,
( h i + h j ) / 2
, as the comparison benchmark. As long as
d i j
is not greater than
( h i + h j ) / 2
are considered to be neighboring particles.
Additionally, an improved auxiliary background grid method is used to accelerate the neighborhood search process. The role of the auxiliary grid is to quickly determine the search scope. The side length of the grid cell is the maximum smoothing length of all the particles; this ensures that all neighboring particle pairs are in the same cell or adjacent cells. During the search process, each particle does not need to match all the particles but only needs to search in its own grid cell or the neighboring grid cells (9 and 27 cells in 2D and 3D, respectively). Thus, the computational complexity of the neighbor search is reduced from
O ( n 2 )
O ( n m )
, where
represent the number of particles in the scene and the average number of particles in the grid cell, respectively. Normally,
m < < n
; thus, the search efficiency is significantly improved.
3.3 Interacting forces
In SPH, the value of physical quantity
at position
is determined by the weighted interpolation sum of its neighboring particles:
A ( x i ) = j m j A j ρ j W x i   - x j , h
m j
represents the mass of particle
ρ j
represents the particle density,
represents the interpolation kernel function , and
represents the smoothing length. Using Equation (5), the pressure on particle
can be calculated as follows:
P ( x i ) = j m j P j ρ j W x i   - x j , h
P j
represents the pressure on particle
. Thus, the net pressure force on particle
is given as:
f i n e t   p r e s s u r e = - P = - j m j P j ρ j W x i   - x j , h = - j V j P j W x i   - x j , h
V j
represents the particle volume. By integrating the pressure over the volume of the particle, the whole pressure force can be determined as follows:
f i p r e s s u r e = - V i P = - V i j V j P j W x i   - x j , h .
This pressure force is not symmetric; thus, we replace
P j
( P i + P j ) / 2
, yielding:
f i p r e s s u r e = - V i j V j P i + P j 2 W x i   - x j , h .
Because we assign different smoothing lengths to the particles, a problem arises: the forces exerted onto adjacent particles differ. For example, the pressure force exerted by particle
onto particle
is given as:
f i j p r e s s u r e = - V i V j P i + P j 2 W ( x i   - x j , h i )
while the pressure force exerted by particle
onto particle
is given as:
f j i p r e s s u r e = - V j V i P i + P j 2 W ( x i   - x j , h j ) .
Because we normally have
h i h j
, the calculation results of the interpolation kernel function in the above equations are different, i.e.,
W ( x i   - x j , h i ) W ( x i   - x j , h j )
. Thus, the forces between the neighboring particles are different, i.e.,
f i j p r e s s u r e f j i p r e s s u r e
, which clearly violates Newton's third law. To solve this problem, we used a combination of their kernel functions. The proposed kernel function
W i j
is computed as the average of the kernel functions of the neighboring particles.
W i j = W ( x i   - x j , ( h i + h j ) / 2 ) .
According to the proposed symmetric kernel function
W i j
, we have
f i j p r e s s u r e = f j i p r e s s u r e
. Similarly, other forces such as the viscous force can be defined as symmetric ones.
4 Experimental results
We tested the proposed ASLSPH algorithm and compared it with state-of-the-art methods, including WCSPH[13], PCISPH[14], IISPH[15], and DFSPH[16,17]. All the experimental results reported in this section were obtained using a desktop computer with an Intel i7-7700k central processing unit and 8 GB of memory.
4.1 Number of neighbors
First, we compared the ASLSPH method and the other four methods with regard to the number of neighbors via a 2D dam-break experiment. As shown in Figure 2, different colors were used to represent the number of neighboring particles for each particle, and a change in color from purple to red indicates that the number of neighboring particles changed from 0 to 20. Compared with that for the other methods, the number of neighboring particles was more evenly distributed for the ASLSPH method; this is because the smoothing length was adaptively adjusted to keep the number of neighbors as consistent as possible.
Additionally, we compared the variance of the number of particle neighbors among the various methods at different timesteps, as shown in Figure 3. The X-axis indicates the timestep (the sampling interval was 20), and the Y-axis indicates the variance of the number of particle neighbors for each method at the corresponding timestep. The variance of the ASLSPH method was significantly smaller than those of the other methods.
4.2 Density field
We also compared the density field of the ASLSPH method and those of the other four methods in the foregoing 2D dam-break experiment. As shown in Figure 4, the density distribution corresponding to the ASLSPH method was significantly more uniform than that of WCSPH. The iterative correction algorithm of PCISPH failed on the left side of the fluid, exhibiting an uneven density distribution. IISPH and DFSPH exhibited uniform density distributions, but both algorithms required substantially more computations than those required by ASLSPH; this is because they needed to solve the complex pressure Poisson equations.
4.3 Visual effects
Next, we compared the visual quality exhibited by the different methods. Various colors were used to represent the smoothing lengths of the particles in the particle view, and a color change from blue to red indicated that the smoothing length changed from small to large.
Figure 5 shows a simulated 3D dam-break experiment with the water column initially located in the middle. The top row presents the results for ASLSPH in the particle view, and the bottom row presents the corresponding surface rendering results. As shown, the simulation effect was highly realistic, with rich details. Figure 6 compares the simulation results of ASLSPH (left), PCISPH (middle), and WCSPH (right). Owing to the lack of neighboring particles on the free surface, the smoothing lengths of the particles on the free surface increased in ASLSPH, and the colors of the particles on the surface gradually changed to green, yellow, and even red. Clearly, more water droplets were captured compared with the case of PCISPH, and even the effect of a rolling wave was captured. WCSPH uses a gas state equation to calculate the pressure, and the fluid settles down quickly, making it difficult to capture the wave curling effect.
As shown in Figure 7, we simulated the effect of a water column falling onto a dragon-shaped obstacle. The simulation results of ASLSPH and IISPH are presented in the top and bottom rows, respectively. Near the free surface, the ASLSPH method effectively adjusted the smoothing length of the particles, improving the calculation accuracy. Compared with IISPH, ASLSPH captured richer surface spray details.
As shown in Figure 8, we simulated the animation result of a bunny-shaped fluid falling into a water tank. The simulation results of ASLSPH and DFSPH are presented in the top and bottom rows, respectively. For the aforementioned reason, the ASLSPH method dynamically adjusted the smoothing length of the particles near the free surface. Initially, DFSPH captured rich fluid details because of its ability to maintain the incompressibility of fluids. However, this capability decreased over time, and DFSPH did not maintain sufficient details (in contrast to ASLSPH), as shown in the last column of Figure 8.
4.4 Simulation efficiency
To evaluate the time efficiencies of the different methods, we simulated the 3D dam-break effect shown in Figure 5 with two different resolutions (particle numbers of 9360 and 89216). Table 1 presents the performance measurements for the first 15s, where Num represents the particle number, △t represents the timestep, Tsim represents the total simulation time, Tns represents the time cost for the neighbor search, and Speed-up represents the acceleration ratio of each method to WCSPH.
Performance measurements of different methods with two simulation resolutions
Model Num t [s] Tsim [s] Tns [s] Speed-up
WCSPH 9360 0.00004 1679.84 486.31 -
PCISPH 0.00130 106.86 29.69 15.72
IISPH 0.00150 82.55 23.27 20.35
DFSPH 0.00150 41.10 11.64 40.87
ASLSPH 0.00150 34.54 11.50 48.63
WCSPH 89216 0.00002 28594.19 8658.68 -
PCISPH 0.00065 1690.96 533.61 16.91
IISPH 0.00100 1233.57 370.13 23.18
DFSPH 0.00150 704.64 209.74 40.58
ASLSPH 0.00150 635.99 231.87 44.96
Although the calculation times of PCISPH, IISPH, and DFSPH in a single timestep were longer than that of WCSPH owing to the iterative calculation of the pressure, larger timesteps could be used (32.5, 37.5, and 37.5, respectively, times that of WCSPH at the low resolution). Consequently, the overall acceleration ratio reached 15.72, 20.35, and 40.87, respectively, in the low-resolution simulation. Owing to the high interpolation accuracy and stability, ASLSPH can also adopt a large timestep. Additionally, similar to WCSPH, ASLSPH does not need to solve the pressure iteratively. ASLSPH achieved an exceedingly high computational efficiency, reaching an acceleration ratio of 48.63 times that of WCSPH. In the high-resolution simulation, the acceleration ratios of PCISPH, IISPH, DFSPH, and ASLSPH were 16.91, 23.18, 40.58, and 44.96 times larger, respectively, than that of WCSPH.
ASLSPH did not perform well with regard to the time required for the neighbor search, even though we used a scheme (see Section 3.2) for acceleration. In the low-resolution mode, the neighbor search times of the five methods accounted for 28.95%, 27.78%, 28.19%, 28.32%, and 33.30% of the total time, respectively (calculated as Tns /Tsim ). The proportion was largest for ASLSPH. In the high-resolution mode, as the number of particles increased, ASLSPH's neighbor search became even more time-consuming, accounting for 36.46% of the total time. This is because the size of the auxiliary grid depends on the maximum smoothing length of the particles, and the efficiency of the neighbor search method in ASLSPH decreases when the smoothing length varies over a large range. Therefore, enhancing the performance of the neighbor search scheme can improve the overall simulation efficiency of ASLSPH.
5 Conclusion and future work
An ASLSPH method is proposed to improve the accuracy of interpolation kernels in SPH by adjusting the smoothing length of particles adaptively. An effective neighboring particle search scheme and a symmetric kernel function are used to increase the time efficiency and maintain the stability of the algorithm, respectively. Thus, a larger timestep can be used in ASLSPH to improve the simulation efficiency. Compared with other state-of-the-art methods, ASLSPH is competitive with regard to both the visual quality and the simulation efficiency.
The ASLSPH method also has shortcomings. When the smoothing length of the particles varies widely, the efficiency of the grid-assisted neighborhood search scheme decreases. In the future, we will investigate better neighbor search schemes to further improve the simulation efficiency. Additionally, we are planning to transplant ASLSPH to a GPU for a large-scale fluid simulation.



Müller M, Solenthaler B, Keiser R, Gross M. Particle-based fluid-fluid interaction. In: Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation-SCA'05. Los Angeles, California, New York, ACM Press, 2005 DOI:10.1145/1073368.1073402


Ihmsen M, Orthmann J, Solenthaler B, Kolb A, Teschner M. SPH fluids in computer graphics. Eurographics, 2014(2): 21–42 DOI:10.2312/egst.20141034


Koschier D, Bender J, Solenthaler B. Smoothed particle hydrodynamics techniques for the physics based simulation of fluids and solids. Eurographics, 2019 DOI:10.2312/egt.20191035


Adams B, Pauly M, Keiser R, Guibas L J. Adaptively sampled particle fluids. ACM Transactions on Graphics, 2007, 26(3): 48 DOI:10.1145/1276377.1276437


Chaniotis A K, Poulikakos D, Koumoutsakos P. Remeshed smoothed particle hydrodynamics for the simulation of viscous and heat conducting flows. Journal of Computational Physics, 2002, 182(1): 67–90 DOI:10.1006/jcph.2002.7152


Lapenta G, Brackbill J U. Control of the number of particles in fluid and MHD particle in cell methods. Computer Physics Communications, 1995, 87(1/2): 139–154 DOI:10.1016/0010-4655(94)00180-a


Miller G, Pearce A. Globular dynamics: a connected particle system for animating viscous fluids. Computers & Graphics, 1989, 13(3): 305–309 DOI:10.1016/0097-8493(89)90078-2


Gingold R A, Monaghan J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 1977, 181(3): 375–389 DOI:10.1093/mnras/181.3.375


Shao S D, Lo E Y M. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface. Advances in Water Resources, 2003, 26(7): 787–800 DOI:10.1016/s0309-1708(03)00030-7


Ren B, Fan H F, Bergel G L, Regueiro R A, Lai X, Li S F. A peridynamics-SPH coupling approach to simulate soil fragmentation induced by shock waves. Computational Mechanics, 2015, 55(2): 287–302 DOI:10.1007/s00466-014-1101-6


Gong K, Shao S D, Liu H, Lin P Z, Gui Q Q. Cylindrical SPH simulations of water entry. Journal of Fluids Engineering, 2018, 141(7): 071303 DOI:10.1115/1.4042369


Kazemi E, Koll K, Tait S, Shao S D. SPH modelling of turbulent open channel flow over and within natural gravel beds with rough interfacial boundaries. Advances in Water Resources, 2020, 140: 103557 DOI:10.1016/j.advwatres.2020.103557


Becker M, Teschner M. Weakly compressible SPH for free surface flows. In: Proceedings of the Acm Siggraph/eurographics Symposium on Computer Animation. New York, Acm Press, 2007 DOI: 10.2312/SCA/SCA07/209-218


Solenthaler B, Pajarola R. Predictive-corrective incompressible SPH. ACM Transactions on Graphics, 2009, 28(3): 1–6 DOI:10.1145/1531326.1531346


Ihmsen M, Cornelis J, Solenthaler B, Horvath C, Teschner M. Implicit incompressible SPH. IEEE Transactions on Visualization and Computer Graphics, 2014, 20(3): 426–435 DOI:10.1109/tvcg.2013.105


Angeles Los, California, York New, Press ACM, 2015, 147–155 DOI:10.1145/2786784.2786796


Bender J, Koschier D. Divergence-free SPH for incompressible and viscous fluids. IEEE Transactions on Visualization and Computer Graphics, 2017, 23(3): 1193–1206 DOI:10.1109/tvcg.2016.2578335


Lastiwka M, Quinlan N, Basa M H. Adaptive particle distribution for smoothed particle hydrodynamics. International Journal for Numerical Methods in Fluids, 2005, 47(10/11): 1403–1409 DOI:10.1002/fld.891


Keiser R, Adams B, Dutré P. Multiresolution particle-based fluids. Swiss Federal Institute of Technology Zurich, Department of Computer Science, 2006 DOI:10.3929/ethz-a-006780981


Angeles Los, California, USA, IEEE, 2008 DOI:10.2312/VG/VG-PBG08/137-146


Orthmann J, Kolb A. Temporal blending for adaptive SPH. Computer Graphics Forum, 2012, 31(8): 2436–2449 DOI:10.1111/j.1467-8659.2012.03186.x


Desbrun M, Gascuel M P. Smoothed particles: a new paradigm for animating highly deformable bodies. Computer Animation and Simulation '96, 1996 DOI:10.1007/978-3-7091-7486-9_5


Ihmsen M, Akinci N, Gissler M. Boundary handling and adaptive time-stepping for PCISPH. The Eurographics Association, 2010 DOI:10.2312/PE/vriphys/vriphys10/079-088


Goswami P, Pajarola R. Time adaptive approximate SPH. In: Proceedings of the 8th Workshop on Virtual Reality Interactions and Physical Simulations, 2011 DOI: 10.2312/PE/vriphys/vriphys11/019-028


Gingold R A, Monaghan J J. Kernel estimates as a basis for general particle methods in hydrodynamics. Journal of Computational Physics, 1982, 46(3): 429–453 DOI:10.1016/0021-9991(82)90025-0


Nelson R P, Papaloizou J C B. Variable smoothing lengths and energy conservation in smoothed particle hydrodynamics. Monthly Notices of the Royal Astronomical Society, 1994, 270(1): 1–20 DOI:10.1093/mnras/270.1.1


Di L, Sigalotti G, Daza J, Donoso A. Modelling free surface flows with smoothed particle hydrodynamics. Condensed Matter Physics, 2006, 9(2): 359 DOI:10.5488/cmp.9.2.359


Olejnik M, Szewc K, Pozorski J. SPH with dynamical smoothing length adjustment based on the local flow kinematics. Journal of Computational Physics, 2017, 348: 23–44 DOI:10.1016/


Ren H L, Zhuang X Y, Rabczuk T. A dual-support smoothed particle hydrodynamics for weakly compressible fluid inspired by the dual-horizon peridynamics. Computer Modeling in Engineering & Sciences, 2019, 121(2): 353–383 DOI:10.32604/cmes.2019.05146


Harada T, Koshizuka S, Kawaguchi Y. Smoothed particle hydrodynamics on GPUs. Computer Graphics International, 2007, 40: 63–70 DOI:10.1142/S0219876207000972


Ihmsen M, Akinci N, Becker M, Teschner M. A parallel SPH implementation on multi-core CPUs. Computer Graphics Forum, 2011, 30(1): 99–112 DOI:10.1111/j.1467-8659.2010.01832.x


Sin F, Bargteil A W, Hodgins J K. A point-based method for animating incompressible flow. In: Proceedings of the 2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation-SCA '09. New Orleans, Louisiana, New York, ACM Press, 2009, 247–255 DOI:10.1145/1599470.1599502


Hernquist L, Katz N. TREESPH-A unification of SPH with the hierarchical tree method. The Astrophysical Journal Letters Supplement Series, 1989, 70: 419 DOI:10.1086/191344


Gong X, Yang J, Zhang S. A parallel SPH method with background grid of adaptive mesh refinement. Chinese Journal of Computational Physics, 2016, 33(2):183-189 DOI:10.3969/j.issn.1001-246X.2016.02.008


Goswami P, Schlegel P, Solenthaler B, Pajarola R. Interactive SPH simulation and rendering on the GPU. In: Proceedings of the 2010 ACM SIGGRAPH/Eurographics Symposium on Computer Animation. 2010, 55–64 DOI:10.2312/SCA/SCA10/055-064


Macklin M, Müller M. Position based fluids. ACM Transactions on Graphics, 2013, 32(4): 1–12 DOI:10.1145/2461912.2461984


Hochstetter H, Kolb A. Evaporation and condensation of SPH-based fluids. In: Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation. Los Angeles, California, New York, NY, USA, ACM, 2017 DOI:10.1145/3099564.3099580


Park S H, Jo Y B, Ahn Y, Choi H Y, Choi T S, Park S S, Yoo H S, Kim J W, Kim E S. Development of multi-GPU–based smoothed particle hydrodynamics code for nuclear thermal hydraulics and safety: potential and challenges. Frontiers in Energy Research, 2020, 8: 86 DOI:10.3389/fenrg.2020.00086


Morikawa D, Senadheera H, Asai M. Explicit incompressible smoothed particle hydrodynamics in a multi-GPU environment for large-scale simulations. Computational Particle Mechanics, 2020, 1–18 DOI:10.1007/s40571-020-00347-0