EUROGRAPHICS 2015 / O. Sorkine-Hornung and M. Wimmer(Guest Editors)

Volume 34 (2015), Number 2

Implicit Formulation for SPH-based Viscous Fluids
Tetsuya Takahashi1,2 Yoshinori Dobashi3,2 Issei Fujishiro4 Tomoyuki Nishita2,5 Ming C. Lin1
1 The University of North Carolina at Chapel Hill, USA
2 UEI Research, Japan

3

Hokkaido University, Japan
4
Keio University, Japan
5
Hiroshima Shudo University, Japan

Figure 1: Viscous fluids simulated with our implicit formulation. Left to right: caramel sauce coiling with a particle view inthe inset; a dragon consisting of particles with different viscosities; melted chocolate buckling with a particle view in the inset.
Abstract
We propose a stable and efficient particle-based method for simulating highly viscous fluids that can generatecoiling and buckling phenomena and handle variable viscosity. In contrast to previous methods that use explicit
integration, our method uses an implicit formulation to improve the robustness of viscosity integration, therefore
enabling use of larger time steps and higher viscosities. We use Smoothed Particle Hydrodynamics to solve the fullform of viscosity, constructing a sparse linear system with a symmetric positive definite matrix, while exploiting thevariational principle that automatically enforces the boundary condition on free surfaces. We also propose a new
method for extracting coefficients of the matrix contributed by second-ring neighbor particles to efficiently solve
the linear system using a conjugate gradient solver. Several examples demonstrate the robustness and efficiencyof our implicit formulation over previous methods and illustrate the versatility of our method.
Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism—Animation

1. Introduction
Smoothed Particle Hydrodynamics (SPH) is becoming increasingly popular for simulating fluids because of its attractive features including automatic conservation of mass, implicit tracking of surfaces with frequent topology changes,
and no need for grid structures or meshes. SPH has beendeveloped in various directions, e.g., enforcing fluid incompressibility [BT07, SP09, ICS∗ 14], handling fluid-fluidand fluid-solid interactions [MSKG05,SSP07,SP08,BTT09,
AIA∗ 12, HLW∗ 12, RLY∗ 14], and improving computationalc 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John
Wiley & Sons Ltd. Published by John Wiley & Sons Ltd.

efficiency and saving memory usage [APKG07, IABT11,
SG11, OK12], and is recognized as a state-of-the-art fluid
solver in computer graphics [IOS∗ 14].
Over the past decades, various SPH methods have been
proposed and used for a variety of fluid effects in the literature. However, most of these SPH methods assume thatfluid is inviscid or slightly viscous; thus an effective SPH
method that can simulate highly viscous fluids has not yet
been established although we see various viscous materials(e.g., honey, caramel sauce, melted chocolate, lava, machin-

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids

ery oils, and bodily fluids) and their characteristic behaviors
on a daily basis. There are two main reasons for this; First,
previous SPH methods, e.g., [MCG03, MSKG05, SSP07],
drop off-diagonal components of viscous stress tensor tosimplify the viscosity term and consequently fail to generate rotational viscous fluid behaviors, such as coiling andbuckling, due to neglect of the boundary condition on freesurfaces, while leading to inaccurate handling of variable
viscosity. This is also true for XSPH [Mon89] and artificial
viscosity [Mon92] as they are essentially same as the Laplacian form in simplifying the viscosity term. Since the simplified term is described by Laplacian operator, in this paper
we refer to this term as Laplacian form of viscosity and theoriginal, unsimplified term as full form of viscosity. Second,previous SPH methods, e.g., [PPLT06,ASP∗ 14], suffer froma restriction on available time steps because of their explicit
viscosity integration schemes. Andrade et al. [ASP∗ 14] pro2

posed a condition on time steps ∆t ≤ 0.1 ρh
8µ (∆t: time step,
ρ: fluid density, h: kernel radius, and µ: dynamic viscosity)
to perform numerically stable simulations of viscous fluidswith the full form, and this condition makes it difficult tosimulate viscous fluids within a reasonable time when higher
viscosity and resolutions (smaller kernel radii) are used.
To address the two problems above, we propose a new
SPH-based method that uses implicit viscosity integration
for the full form for robustly simulating highly viscous fluids. Our method offers the following advantages:• It is efficient, allowing for use of larger time steps andfiner spatial resolutions than explicit integration.
• It is robust and stable, even with large time steps and highviscosities.
• It can generate coiling and buckling phenomena and handle variable viscosity.
We exploit the variational principle that automaticallyenforces the boundary condition on free surfaces to derive
our implicit formulation, constructing a sparse linear system with a symmetric positive definite matrix. To efficiently
solve the linear system, we also propose a novel method
for extracting coefficients of the matrix that includes contributions from first-ring neighbor particles and second-ringneighbor particles (neighbor particles’ neighbor particles).
Figure 1 demonstrates viscous fluids, simulated using our
implicit formulation.
In Eulerian methods, we can easily discretize the fullform (divergence of Jacobian of velocity) at a time using
finite difference due to the staggered stress arrangement as
in [BB08]. On the other hand, in SPH, we first need to compute Jacobian of velocity, and then compute divergence ofJacobian of velocity, separately applying SPH formulations
to both steps. Although computing these two steps is easy
for explicit integration [PPLT06, ASP∗ 14], significant complexity is involved with implicit integration, and we need totake into account contributions from both of first-ring andsecond-ring neighbors to construct a linear system. Becauseof this complexity unique to SPH discretization and derivation of the implicit formulation, no work has been proposed
regardless of the apparent simplicity of the idea in adopting
implicit integration and demand for robust SPH-based simulators. To the best of our knowledge, our method is the first
SPH method that uses implicit integration for the full formof viscosity, and is also the first method that extracts matrix
coefficients contributed by second-ring neighbors.
2. Related Work
Eulerian viscous fluids. Carlson et al. [CMVHT02] first enabled stable simulations of highly viscous fluids with freesurfaces by solving the Laplacian form of viscosity using
implicit integration. Later, Rasmussen et al. [REN∗ 04] proposed an implicit-explicit scheme for the full form of viscosity to correctly handle variable viscosity at the expenseof numerical stability. Batty and Bridson [BB08] proposed afully implicit viscosity integration scheme for the full form,
making it possible to take larger time steps, handle variable
viscosity, and generate coiling and buckling. This method
was extended by Batty and Houston [BH11] for an adaptive tetrahedral fluid simulator. In the work of Stomakhin etal. [SSJ∗ 14], a viscosity term was also solved using implicit
integration in the framework of Material Point Method.
Lagrangian viscous fluids. We categorize Lagrangian
methods into five groups: Lagrangian Finite Element Methods (Lagrangian FEM), dimensionally reduced discrete
methods, spring-based methods, deformation-based methods, and SPH methods. Our method belongs to the SPH
methods.
Lagrangian FEM has been used to accurately simulate viscous fluids, and various developments have been done in theliterature. Bargteil et al. [BWHT07] proposed an efficient
remeshing method to reduce the cost of time-consuming
remeshing process, and Wojtan and Turk [WT08] improved
the remeshing method of Bargteil et al. [BWHT07]. Wicke
et al. [WRK∗ 10] proposed a local remeshing method to keep
the number of tetrahedra small. Clausen et al. [CWSO13]
proposed a Lagrangian FEM that can handle elastic, plastic,
and fluid materials in a unified manner.
To accurately simulate viscous threads and sheets, Bergou
et al. [BAV∗ 10] and Batty et al. [BUAG12] proposed dimensionally reduced discrete methods and generated coiling andbuckling, respectively, limiting the dimension of materials
that they can simulate.
Spring-based methods have been used because of its
conceptual simplicity. Miller and Pearce [MP89] and Terzopoulos et al. [TPF91] proposed a spring-based model thatcomputes repulsion and attraction forces between particles.
Clavet et al. [CBP05] extended this model to simulate materials that exhibit elasticity, viscosity, and plasticity. Takahashi et al. [TNF14] also simulated such materials in a unified framework of Position-based dynamics.
c 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids

Gerszewski et al. [GBB09] proposed a deformation-based
method that approximates motions of neighbor particles
based on deformations of particle configurations for reproducing elastoplastic materials. Their method was extended
by Zhou et al. [ZLKW13] to improve its numerical stability using implicit integration and by Jones et al. [JWJ∗ 14] tohandle varying mass materials.
Desbrun and Gascuel [DG96] used SPH to simulate viscous materials. Müller et al. [MCG03] proposed the Laplacian form of viscosity and simulated slightly viscous fluid.
The Laplacian form was also used in [MSKG05, SSP07].
Paiva et al. [PPLT06] proposed the full form of viscosityand accurately simulated viscous fluids. The full form wasalso used in the work of Andrade et al. [ASP∗ 14]. Rafiee etal. [Raf07] presented a method based on a Maxwell model
to simulate coiling of viscoelastic fluids. Dagenais et al.
[DGP12] simulated viscous fluid motions by adding extra
forces that move particles to their original positions. In Astrophysics, Monaghan [Mon97] and Laibe and Price [LP12]
used pair-wise implicit formulations for artificial viscosity.
Their methods consider only pair-wise particles for implicit
integration, and thus can be numerically unstable. In Computational Mechanics, Fan et al. [FTZ10] proposed an implicit scheme for simulating viscous fluids using SPH. However, they computed viscous stress with gradient of velocity,
not Jacobian of velocity, and hence their method is essentially equivalent to a method that uses the Laplacian form.
Paiva et al. [PPLT06] and Andrade et al. [ASP∗ 14] solved
the full form of viscosity. However, our method differs fromtheirs in that theirs used explicit integration whereas ours
uses implicit integration to improve the robustness of thesimulation and enable use of larger time steps with higher
viscosities and finer spatial resolutions than these methods.
3. Fundamentals for Simulating Viscous Fluids
Formulations. We aim to simulate incompressible, highly
viscous fluids using SPH. In the Lagrangian setting, theNavier-Stokes equations for particle i can be described as
ρi

dui
ρ
= −∇pi + ∇ · si + i Fext
,
dt
m i


si = µi ∇ui + (∇ui )T ,

(1)
(2)

where ρi denotes density of particle i, t time, ui = [ui , vi , wi ]T
velocity, pi pressure, si viscous stress tensor, m mass (weuse a constant mass for all particles), Fext
external force,iand µi dynamic viscosity. The combination of the second
term on the right side in Eq. (1), and Eq. (2) is the fullform of viscosity and is required to handle variable viscosity [REN∗ 04] and rotational behaviors [BB08]. We separate the terms in Eq. (1) to independently solve them, taking
the standard approach of operator splitting as with Eulerian
methods [BB08, BH11], and enforce fluid incompressibility using a particle-based incompressible fluid solver, e.g.,[SP09, ICS∗ 14].
c 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.

To generate rotational behaviors by solving Eq. (1), we
need to consider the boundary condition that there is no traction on free surfaces [BB08]. In other words, as a boundary condition for Eq. (1), we must satisfy (−pi I + si )ni = 0(ni : normal to the free surface). As with [BB08,SSJ∗ 14], we
decouple these terms and independently enforce the boundary condition for pressure −pi Ini = 0 (we enforce this bysetting pressuresto 0), and for viscos of surface particlesT
ity si ni = µi ∇ui + (∇ui ) ni = 0. Since our method isbased on the variational principle, the boundary condition
for viscosity is automatically enforced by solving Eq. (1)(see e.g., [BB08, SSJ∗ 14]).
Algorithm. First, we apply only external force Fext
i andobtain first intermediate velocity u∗i . Then, we solve viscosity using our implicit formulation (see § 4) and obtain second
intermediate velocity u∗∗
i with intermediate viscous stresstensor s∗∗
:
i
∆t
∗
u∗∗
∇ · s∗∗
i = ui +
i ,
ρi


∗∗
∗∗ T
s∗∗
.
i = µi ∇ui + ∇ui

(3)
(4)

Next, to enforce fluid incompressibility, we compute pressure pi using a particle-based fluid solver with the boundp
ary condition for pressure, and compute pressure force Fi[Mon92]. Finally, particle velocity ui and position xi =
[xi , yi , zi ]T are integrated using Euler-Cromer scheme. We
summarize a full procedure of our method in Algorithm 1.
Algorithm 1 Procedure of our method
1: // j: neighbor particle of i
2: // Wi j : kernel with a kernel radius h
3: for all particle i do
4:find neighbor particles
5: for all particle i do
6:
apply external force u∗i = uti + ∆tFext
i /m
7: for all particle i do
8:solve viscosity using Eqs. (3) and (4) // § 4
9: for all particle i do
10:compute pi using a particle-based fluid solver
11: for all particle i do


12:

p

compute Fi = −m2 ∑ j

p
pi
+ ρ2j
ρ2i
j

∇Wi j

13: for all particle i dop
14:
integrate particle velocity ut+1
= u∗∗
i + ∆tFi /mi15:

integrate particle position xt+1
= xti + ∆tut+1ii

Unlike [BB08] that enforces fluid incompressibility twice,
we enforce fluid incompressibility once in one simulation
step as in [BH11] because our method can sufficiently maintain incompressibility over all simulation steps (when fluid
incompressibility was enforced with a tolerable density error of 0.1%, the maximum density deviation was lower than0.5%), and we observed indiscernible differences between
enforcing incompressibility once and twice.

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids

4. Implicit Formulation for Full Form of Viscosity
In this section, we first describe how to solve the viscosity term using implicit integration in the SPH framework
while constructing a linear system (§ 4.1). Next, we explain
sparsity of the coefficient matrix (§ 4.2), and our solver andcoefficient extraction method for the system (§ 4.3). Then,we give implementation details and show our algorithm forsolving our implicit formulation (§ 4.4).
4.1. Implicit Integration for Full Form of Viscosity
We derive an SPH-based implicit viscosity formulation from
Eqs. (3) and (4). Hereafter, we drop the symbol ∗∗ for readability.
Since the boundary condition for viscosity is automatically enforced because of the variational principle [BB08,
SSJ∗ 14], we directly discretize Eqs. (3) and (4) using implicit integration in the SPH framework:
!
sj
si
∗
ui = ui + m∆t ∑+ 2 ∇Wi j , (5)
ρ2i
ρj
j


si = µi ∑ V j (u j − ui )∇WiTj + ∇Wi j (u j − ui )T . (6)jBy substituting si in Eq. (6) into Eq. (5) and arranging theterms in these equations, we obtain an implicit formulation:
ui + m̂ ∑ µ̂i Qi j + µ̂ j Q jk ∇Wi j = u∗i , (7)jqi j,xy
qi j,xz
 , (8)2 ∑ j ai j,y vi j
qi j,yz
qi j,yz
2 ∑ j ai j,z wi j
qi j,xy = ∑ ai j,y ui j + ai j,x vi j , qi j,xz = ∑ ai j,z ui j + ai j,x wi j ,

2 ∑ j ai j,x ui j
qi j,xy
Qi j = qi j,xzjj
qi j,yz = ∑ ai j,z vi j + ai j,y wi j ,jwhere m̂ = m∆t, µ̂i = µi /ρ2i , k is a neighbor particle of j, ai j = [ai j,x , ai j,y , ai j,z ]T = V j ∇Wi j =V j [∇Wi j,x , ∇Wi j,y , ∇Wi j,z ]T , ui j = ui − u j , vi j = vi − v j , andwi j = wi − w j . Qi j is symmetrical due to the symmetricalproperty of si . This implicit formulation Eq. (7) is a linear
system and can be rewritten in a matrix form as CU = U∗
(C is a coefficient matrix and U = [. . . , ui , vi , wi , . . .]T ). Let
N denote the number of fluid particles, the size of C is
3N × 3N, and that of U is 3N × 1. We assign serial numbers(id) to particles to specify locations in a coefficient matrix.
4.2. Sparsity of Coefficient Matrix
In SPH, particle i interacts with its neighbor particles within
its kernel radius, namely inside of its kernel sphere. Since
our formulation requires viscous stress of particles i and j
to compute velocity update for particle i using Eq. (5), notonly i’s neighbors Ji (a set of particle j) but j’s neighbors Ki j

Mi

Mi

h

h
j∗∗

xi

2h
j∗∗

j∗

k∗(i)

i
Ki j∗∗Ji

Ki j∗

Ki j∗∗

2h

xi

Ji

j∗
k∗∗
Ki j∗

Figure 2: Illustration for first- and second-ring neighbors.
Left: particle i (green) directly interacts with i’s first-ringneighbors Ji (particle j colored in orange) inside of i’s kernel sphere with its radius h shown as a green circle. Right:
particle j∗ (purple) has neighbor particles (within j∗ ’s kernel sphere shown as a purple circle), which are second-ringneighbors Ki j∗ for particle i. Possible second-ring neighbors
for particle i exist within the largest Minkowski sum (Mi ) ofJi and Ki j , which is a sphere shown as a red circle, and Mi ’s
center and radius are xi and 2h, respectively.
(a set of particle k) must be taken into account. In short, weneed to include contributions from first-ring neighbors Ji andsecond-ring neighbors Ki j to compute i’s velocity update andconstruct a linear system. This setup is illustrated in Figure
2. Assuming that particles i and j are spherically surroundedby others, i and j generally have 30-40 first-ring neighbors
within their kernel radius h [SP09]. Since particle j can exist anywhere within i’s kernel sphere, the total number of
second-ring neighbors of particle i without overlaps can be240-320. This is because the maximum (total) number of
second-ring neighbors of i without overlaps is smaller than
the number of particles within the largest Minkowski sumof i’s kernel sphere and js’ kernel spheres, and the largest
Minkowski sum Mi is a sphere whose center is xi and radius
is 2h. Since all particles in Ji and Ki j are included in Mi ,each particle can interact with up to 320 particles. Hence,
non-zero values for each velocity component can be 960, as
there are 320 particles in Mi and each has three velocity components. Although this number is much larger than the number of non-zero values in grid-based methods, and particlebased methods that involve only first-ring neighbors, our coefficient matrix is still a sparse matrix.
4.3. Solver and Coefficient Extraction
Solver. Since the linear system constructed by our formulation is sparse and also symmetric positive definite, we use
a conjugate gradient (CG) solver (though we also tested Jacobi method and Modified Incomplete Cholesky Conjugate
Gradient (MICCG), they did not work well. See § 6). Unlike
the CG method described in [Ihm13], which repeatedly computes matrix coefficients in the CG algorithm by performing
extra particle scans without storing the coefficients, we explicitly construct and preserve a coefficient matrix, namely
extract all coefficients in the matrix. This approach enables
us to efficiently perform the CG method without extra loops,
c 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids


T 
uk
cuk ui∑  cvk ui   vk  = u∗i ,kwk
cwk ui
cui ui = 1 + m̂µ̂i 2ωi j,x αi j,x + ωi j,y αi j,y + ωi j,z αi j,z ,

cvi ui = m̂µ̂i ωi j,y αi j,x ,cwi ui = m̂µ̂i ωi j,z αi j,x ,
cu j ui = m̂ −µ̂i (2ai j,x ωi j,x + ai j,y ωi j,y + ai j,z ωi j,z ) +Figure 3: Numerical stability test with different combinations of time steps and viscosities. (a) initial state, (b)explicit integration [ASP∗ 14] with ∆t = 5.0 × 10−6 s andµ = 1, 000.0 kg/(s·m), (c) explicit integration [ASP∗ 14] with∆t = 1.3 × 10−3 s and µ = 1, 000.0 kg/(s·m), (d) implicit
integration (our method) with ∆t = 1.3 × 10−3 s and µ =1, 000.0 kg/(s·m), (e) explicit integration [ASP∗ 14] with∆t = 5.0 × 10−6 s and µ = 50, 000.0 kg/(s·m), and (f) implicit integration (our method) with ∆t = 1.0 × 10−4 s andµ = 50, 000.0 kg/(s·m).

µ̂ j (2∇Wi j,x α jk,x + ∇Wi j,y α jk,y + ∇Wi j,z α jk,z ) ,
cv j ui = m̂ −µ̂i ai j,x ωi j,y + µ̂ j ∇Wi j,y α jk,x ,
cw j ui = m̂ −µ̂i ai j,x ωi j,z + µ̂ j ∇Wi j,z α jk,x ,cuk ui = −m̂ ∑ µ̂ j (2∇Wi j,x a jk,x + ∇Wi j,y a jk,y + ∇Wi j,z a jk,z ),j(10)
cvk ui = −m̂ ∑ µ̂ j ∇Wi j,y a jk,x ,(11)

cwk ui = −m̂ ∑ µ̂ j ∇Wi j,z a jk,x ,

(12)

j

j

thereby improving the performance of the solver. In addition to this advantage, extracting all coefficients allows us
to use preconditioning techniques, e.g., algebraic multigrid(AMG), and external libraries, separate solver code fromothers to improve programming maintainability, and takefull advantage of GPGPU, parallelizing matrix-vector multiplications (in both rows and columns).
After a coefficient matrix is constructed, we solve the linear system with the CG method. We terminate iterations in
the CG algorithm when a relative residual becomes smaller
than a convergence criterion η.
Coefficient extraction. By substituting Qi j in Eq. (8), we
can rewrite Eq. (7) for x component of ui , ui as
ui + m̂ ∑µ̂i 2∇Wi j,x ∑ ai j,x ui j +jj∇Wi j,y ∑(ai j,y ui j + ai j,x vi j ) + ∇Wi j,z ∑(ai j,z ui j + ai j,x wi j )jj
+µ̂ j 2∇Wi j,x ∑ a jk,x u jk + ∇Wi j,y ∑(a jk,y u jk + a jk,x v jk )where αi j = [αi j,x , αi j,y , αi j,z ]T = ∑ j ai j and ωi j =[ωi j,x , ωi j,y , ωi j,z ]T = ∑ j ∇Wi j . We use cui ui to denote a coefficient of ui to ui , and cvi ui a coefficient of vi to ui , andsimilarly define other coefficients. The other components (yand z) and 2D version can be straightforwardly derived. Byscanning particles i, j, and k, we extract all coefficients forthe linear system (see Appendix A for details).
4.4. Implementation Details and Algorithm
We summarize procedures of our implicit formulation in Algorithm 2. To handle solid boundaries, we use solid particles arranged on object surfaces. When fluid particles collide with solid particles, we use explicit viscosity integration for fluid particles with low viscosity while using Dirichlet boundary condition similar to [BB08, SSJ∗ 14, ASP∗ 14],
namely setting averaged solid particle velocities usolid to
fluid particles if viscosity of the fluid particles is higher thana criterion µDirichlet . As a structure of a sparse matrix, we usecompressed sparse row (CSR) and reserve sufficient memory for each particle so that the matrix construction can beparallelized over particle i.

Algorithm 2 Algorithm for solving viscosity
1: assemble the matrix // see Appendix A
2: solve the linear system with CG+∇Wi j,z ∑(a jk,z u jk + a jk,x w jk )  = u∗i . (9) 3: for all fluid particle i dok4:
if µDirichlet < µi ∧ neighbor solid particle exists then5:
enforce solid boundary condition ui = usolidThen, we further convert Eq. (9) into the following equation to straightforwardly extract coefficients
cui ui , cvi ui , cwi ui , cu j ui , cv j ui , cw j ui , cuk ui , cvk ui , and cwk ui :5. Results

T 



T 
Implementation. We implemented our method in C++ andcui uiui
cu j uiuj cvi ui   vi  + ∑  cv j ui   v j  +parallelized it using Open MP 2.0. We adopted SPH kerj
nels proposed in [MCG03] and used surface tension force
cwi uiwi
cw j uiwjkkc 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids

Table 1: Simulation parameters and performance. N: the number of particles, µ (kg/(s · m)): dynamic viscosity of particles,∆t (s): time step, and t visc (s) and t total (s): simulation time for viscosity and total simulation time per frame, respectively.
Viscosity form

3 (b)
3 (c)
3 (d)
3 (e)
3 (f)
5
6 (a), (c)6 (b), (d)7 (a), (c)7 (b), (d)

Full
Full
Full
Full
Full
Full
Laplacian
Full
Full
Full

Integration
Explicit
Explicit
Implicit
Explicit
Implicit
Implicit
Implicit
Implicit
Implicit
Implicit

N
5.5k
5.5k
5.5k
5.5k
5.5k
47.4k
up to 19.6k
up to 46.5k
up to 34.6k
up to 23.4k

presented in [BT07]. We adopted Implicit Incompressible
SPH (IISPH) [ICS∗ 14] as an incompressible fluid solver andused the fluid solid coupling method proposed in [AIA∗ 12].
We used a variant of the z-index neighbor search method
presented in [IABT11]. In addition to our viscosity formulation, we also used artificial viscosity to stabilize simulations [Mon92].

µ∆t1,000.0
1,000.0
1,000.0
50,000.0
50,000.0
up to 800.0
600.0
600.0
100.0
600.0

−6

5.0 × 10
1.3 × 10−3
1.3 × 10−3
5.0 × 10−6
1.0 × 10−4
1.0 × 10−4
2.0 × 10−4
2.0 × 10−4
2.0 × 10−4
2.0 × 10−4

t visct total

6.3
0.0
15.0
2.9
361.4
815.2
20.4
287.1
135.4
116.3

51.7
0.5
15.4
40.1
366.1
836.1
59.3
314.1
142.0
120.7

60
50
40
Time (s)FigureMethod of Andrade et al.
30

Our method

20
10
0
0

50

100

150

200

FramesSetting. We executed all the scenes on a PC with a 4core Intel Core i7 3.40 GHz CPU and RAM 16.0 GB, andrendered all the figures using a physically-based renderer,
Mitsuba. We used fixed time steps and set a convergence criterion as η = 1.0 × 10−4 . Simulation parameters and performance are listed in Table 1, where the surface reconstruction
and rendering are not included in performance measurement.
5.1. Numerical Stability
To verify the numerical stability of our implicit formulation
over the previous method [ASP∗ 14] that uses explicit viscosity integration, we performed a simple test, as shown in Figure 3, where a viscous bunny was dropped onto the ground,and the initial state is shown in Figure 3 (a). In this scene,2

we chose ∆t satisfying the viscosity condition (∆t ≤ 0.1 ρh
8µ )given in [ASP∗ 14], and used IISPH as a fluid solver of[ASP∗ 14] for fair comparison. The method of [ASP∗ 14] cangenerate a plausible behavior of the bunny with a small time
step and low viscosity, as shown in Figure 3 (b). However,when a large time step or high viscosity is used, their method
easily fails to simulate the bunny, as shown in Figures 3 (c)and (e). By contrast, our implicit method successfully simulates the bunny with a much larger time step (Figure 3 (d)),
and with a large time step and high viscosity (Figure 3 (f)).

Figure 4: Performance profile for Figures 3 (b) and (d).

Figure 5: A dragon consisting of particles with different viscosities from 0.0 to 800.0 kg/(s·m). Light (dark) green particles represent low (high) viscosity.
method of [ASP∗ 14], and our computation is 3.4 times fasterthan theirs although per step cost of our implicit method ismore expensive than that of explicit methods. Since time
steps for [ASP∗ 14] are restricted by the viscosity condition above, our method can be more advantageous for scenes
with higher viscosities and finer resolutions.
5.3. Variable Viscosity
Figure 5 illustrates an example of a dragon consisting ofparticles with different viscosities from 0.0 (light green) to
800.0 kg/(m·s) (dark green). Particles flow at different rates
depending on particle velocities.

5.2. Performance
We compared performance of our method and the previous
method [ASP∗ 14] in Figure 4 using the bunny scene shownin Figures 3 (b) and (d). Because of our robust implicit formulation, we can take a 260.0 times larger time step than the

5.4. Buckling and Coiling
We performed a buckling test to demonstrate that our implicit formulation can generate buckling of a viscous material while the Laplacian form with implicit integration failsc 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids

6. Discussions and Limitations
Robustness. Our implicit formulation significantly improves the robustness of viscosity integration and allows
us to avoid the time step restriction for explicit integration [ASP∗ 14]. Our method introduces a much more relaxed
restriction on time steps, which depends only slightly on
viscosity and spatial resolutions. Consequently, our method
may not generate plausible fluid behaviors, when very large
time steps and very high viscosity and resolutions are used,even though an implicit formulation for the Laplacian form
(which is unconditionally stable) can perform stable simulations with the same condition.

Figure 6: Buckling test for comparison of our method with
the Laplacian form using implicit integration. (a) Laplacian
form with meshes. (b) our method with meshes. (c) Laplacian
form with particles. (d) our method with particles.

Figure 7: Coiling test for different viscosities. (a) Low viscosity with meshes. (b) High viscosity with meshes. (c) Low
viscosity with particles. (d) High viscosity with particles.
Light (dark) green particles represent low (high) viscosity.

to generate a buckling phenomenon. Figure 6 illustrates a
dropped chocolate onto the ground, where (a) and (c) are
simulated with the Laplacian form, and (b) and (d) are with
our formulation. The chocolate simulated with the Laplacian
form does not bend, and exhibits unnatural fluid flows regardless of its high viscosity. By contrast, our implicit formulation successfully generates natural chocolate buckling.
We also performed a coiling test with different viscosities in Figure 7. In this scene, we poured caramel sauce withlow and high viscosity onto the ground. While caramel saucewith low viscosity ((a) and (c)) behaves similar to the results generated with the Laplacian form in Figures 6 (a) and(c), caramel sauce with high viscosity exhibits a coiling phenomenon ((b) and (d)).
c 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.

From our experiments, we deduce that one decisive factor for this weaker robustness is due to second-ring neighbors. To examine the factor, we consider solving differential equations with Laplacian operator by using SPH formulations for Laplacian, and also solving the equations byseparating Laplacian into divergence and gradient, applying SPH formulations for both of divergence and gradient operators, and including second-ring neighbors, because
β∇2 φ ≈ ∇·(β∇φ) (β and φ are arbitrary quantities), and thisdecomposition is numerically similar to the relation of the
Laplacian form µ∇2 u and full form ∇ · (µ∇u + µ(∇u)T ).
We actually tested these with the Laplacian form of viscosityand heat equations, and as expected, we observed numericalinstability with larger time steps and higher coefficients andresolutions when Laplacian operator was decomposed while
we were able to perform stable simulations with undecomposed Laplacian operator under the same condition.
Solver. In addition to CG, we tested Jacobi method and
MICCG to solve our linear system. Although we were able
to solve the system using the two solvers, there are a few
problems to note. First, Jacobi method worked for slightly
viscous fluids with smaller time steps under low spatial
resolutions. However, since the convergence rate of Jacobi
method is slow, Jacobi method failed to solve the systemwith large time steps and high viscosity and spatial resolutions that can make a coefficient matrix ill-conditioned. In
contrast to Jacobi method, MICCG worked for such an illconditioned matrix, and the rate of convergence was actually faster than our CG solver. However, MICCG requires
Cholesky factorization that is costly (and difficult to parallelize), especially with our matrix that includes a great number of non-zero values (see § 4) although serial forward andback substitution steps are not so slow compared to othersteps in the CG algorithm. In our experiments, Cholesky factorization occupied more than 90% of computational time
for our viscosity solver including coefficient extractions,and therefore our non-preconditioned CG solver was morethan 6 times faster than MICCG. A similar behavior thatthe incomplete Cholesky preconditioner performed worse
due to the cost of constructing a preconditioning matrix forparticle-based methods (many non-zeros) was also reported
in [HLW∗ 12].

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids

Performance. Solving our viscosity formulation generally occupies more than 90% of the whole computationaltime, and this weakens one of our advantages of efficientcomputation with larger time steps. Our viscosity formulation consists of extracting coefficients through triple loops
and solving a linear system, and they occupy around 30%
and 70% of computational time, respectively, for Figure 3(f). To accelerate the speed of coefficient extractions, an efficient algorithm for extracting coefficients involving secondring neighbors would be useful. One possible fast algorithm
is to avoid triple loops when extracting coefficients, andwe can avoid them by using precomputations for particle
k. However, these computations require additional storage
and scans over the storage, making coefficient extractions
more complex. Consequently, benefits of the precomputations over our method can be lost, or extracting coefficientswith the precomputations can be more costly than ours. As
for solving a linear system, using a faster solver or a low costand effective preconditioner would be helpful.
Memory. Preserving a coefficient matrix requires a large
memory (e.g., 12 GB memory for 500k particles, due to 1kof 8 byte double values for 3 velocity components of 500k
particles). However, this is not a big issue with current memory capacity, given advantages explained in § 4.3.

GPGPU techniques to accelerate constructing and solving
the linear system. In particular, our explicit matrix preservation allows us to take full advantage of dynamic parallelism,
fully parallelizing matrix-vector multiplications. Additionally, finding better solvers, efficient and effective preconditioners as well as fast coefficient extraction methods would
be promising. Similar to improving memory and computational efficiency for enforcing incompressibility in SPH fluids, using multi-sized particles, domain decomposition, andbackground grids could also further optimize performance.
Our coefficient extraction method can be applied to not
only SPH, but also other point-based methods that involve
second-ring neighbors (IISPH [ICS∗ 14] and mesh smoothing [DMSB99]). Extending our method to accelerate suchproblems could also lead to interesting future research directions.
Acknowledgements
This work is supported in part by JASSO for Study Abroad,
JST CREST, JSPS KAKENHI (Grant-in-Aid for Scientific
Research(A)26240015), U.S. National Science Foundation,and UNC Arts and Sciences Foundation. We would liketo thank Ryoichi Ando and anonymous reviewers for theirvaluable suggestions and comments.

Scalability. Another issue to note is that the size of a matrix grows proportionally to the number of particles. It has ascaling factor larger than grid-based methods and particlebased methods that involve only first-ring neighbors. On
larger scenes with up to 200k particles, the memory usage
and computational cost for coefficient extraction increased
at nearly 1.0x linear-scale with respect to the number ofparticles – similar to grid-based methods and particle-based
methods involving only first-ring neighbors – whereas computational cost for solving linear system increased at around1.1x scale. This is because increased number of particles
requires more CG iterations and thus more computational
time. Although computational cost might increase superlinearly when more particles are used, using AMG preconditioners (which cannot be used without explicit matrix preservation) can potentially address this problem.

References

7. Conclusion and Future Work

[BH11] BATTY C., H OUSTON B.: A simple finite volume method
for adaptive viscous liquids. In Proceedings of the 2011 ACM
SIGGRAPH/Eurographics Symposium on Computer Animation(2011), pp. 111–118. 2, 3

We proposed a new SPH-based implicit formulation for thefull form of viscosity. Our method enables efficient andstable viscous fluid simulations with larger time steps andhigher viscosities and resolutions than previous methods thatuse explicit integration while handling variable viscosity andgenerating coiling and buckling. We additionally presented anovel coefficient extraction method for a sparse matrix thatinvolves second-ring neighbors to efficiently solve a linear
system with a CG solver. By taking advantage of our implicit
formulation and coefficient extraction method, we achieved
an accelerated performance by a factor of 3.4.
For future work, we plan to implement our method using[AIA∗ 12] A KINCI N., I HMSEN M., A KINCI G., S OLENTHALER
B., T ESCHNER M.: Versatile rigid-fluid coupling for incompressible SPH. ACM Transactions on Graphics 31, 4 (2012),
62:1–62:8. 1, 6
[APKG07] A DAMS B., PAULY M., K EISER R., G UIBAS L. J.:
Adaptively sampled particle fluids. ACM Transactions on Graphics 26, 3 (2007). 1[ASP∗ 14] A NDRADE L UIZ F. D . S., S ANDIM M., P ETRONETTO
F., PAGLIOSA P., PAIVA A.: SPH fluids for viscous jet buckling.
In Proceedings of SIBGRAPI 2014 (2014), pp. 65–72. 2, 3, 5, 6,
7[BAV∗ 10] B ERGOU M., AUDOLY B., VOUGA E., WARDETZKY
M., G RINSPUN E.: Discrete viscous threads. ACM Transactions
on Graphics 29, 4 (2010), 116:1–116:10. 2[BB08] BATTY C., B RIDSON R.: Accurate viscous free surfaces
for buckling, coiling, and rotating liquids. In Proceedings of the2008 ACM SIGGRAPH/Eurographics Symposium on Computer
Animation (2008), pp. 219–228. 2, 3, 4, 5

[BT07] B ECKER M., T ESCHNER M.: Weakly compressible
SPH for free surface flows. In Proceedings of the 2007 ACM
SIGGRAPH/Eurographics Symposium on Computer Animation(2007), pp. 209–217. 1, 6
[BTT09] B ECKER M., T ESSENDORF H., T ESCHNER M.: Direct
forcing for lagrangian rigid-fluid coupling. IEEE Transactions
on Visualization and Computer Graphics 15, 3 (2009), 493–503.
1[BUAG12] BATTY C., U RIBE A., AUDOLY B., G RINSPUN E.:
Discrete viscous sheets. ACM Transactions on Graphics 31, 4(2012), 113:1–113:7. 2

c 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids[BWHT07] BARGTEIL A. W., W OJTAN C., H ODGINS J. K.,
T URK G.: A finite element method for animating large viscoplastic flow. ACM Transactions on Graphics 26, 3 (2007). 2[CBP05] C LAVET S., B EAUDOIN P., P OULIN P.: Particle-based
viscoelastic fluid simulation. In Proceedings of the 2005 ACM
SIGGRAPH/Eurographics Symposium on Computer Animation(2005), pp. 219–228. 2[CMVHT02] C ARLSON M., M UCHA P. J., VAN H ORN III
R. B., T URK G.: Melting and flowing. In Proceedings of the2002 ACM SIGGRAPH/Eurographics Symposium on Computer
Animation (2002), pp. 167–174. 2[CWSO13] C LAUSEN P., W ICKE M., S HEWCHUK J. R.,
O’B RIEN J. F.: Simulating liquids and solid-liquid interactions
with lagrangian meshes. ACM Transactions on Graphics 32, 2(2013), 17:1–17:15. 2[DG96] D ESBRUN M., G ASCUEL M.-P.: Smoothed particles: A
new paradigm for animating highly deformable bodies. In Proceedings of the Eurographics Workshop on Computer Animationand Simulation (1996), pp. 61–76. 3[DGP12] DAGENAIS F., G AGNON J., PAQUETTE E.:
A
prediction-correction approach for stable SPH fluid simulationfrom liquid to rigid. In Proceedings of the Computer Graphics
International 2012 (2012). 3[DMSB99] D ESBRUN M., M EYER M., S CHRÖDER P., BARR
A. H.: Implicit fairing of irregular meshes using diffusion andcurvature flow. In Proceedings of the 26th Annual Conference
on Computer Graphics and Interactive Techniques (1999), SIGGRAPH ’99, pp. 317–324. 8[FTZ10] FAN X.-J., TANNER R., Z HENG R.: Smoothed particle hydrodynamics simulation of non-Newtonian moulding flow.
Journal of Non-Newtonian Fluid Mechanics 165, 5-6 (2010),
219–226. 3[GBB09] G ERSZEWSKI D., B HATTACHARYA H., BARGTEIL
A. W.: A point-based method for animating elastoplastic solids.
In Proceedings of the 2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (2009), pp. 133–138. 3[HLW∗ 12] H E X., L IU N., WANG G., Z HANG F., L I S., S HAO
S., WANG H.: Staggered meshless solid-fluid coupling. ACM
Transactions on Graphics 31, 6 (2012), 149:1–149:12. 1, 7
[IABT11] I HMSEN M., A KINCI N., B ECKER M., T ESCHNER
M.: A parallel SPH implementation on multi-core CPUs. Computer Graphics Forum 30, 1 (2011), 99–112. 1, 6
[ICS∗ 14] I HMSEN M., C ORNELIS J., S OLENTHALER B., H OR VATH C., T ESCHNER M.: Implicit incompressible SPH. IEEE
Transactions on Visualization and Computer Graphics 20, 3(2014), 426–435. 6, 8
[Ihm13] I HMSEN M.: Particle-based Simulation of Large Bodies
of Water with Bubbles, Spray and Foam. PhD thesis, University
of Freiburg, 2013. 4[IOS∗ 14] I HMSEN M., O RTHMANN J., S OLENTHALER B.,
KOLB A., T ESCHNER M.: SPH fluids in computer graphics. In
EUROGRAPHICS 2014 State of the Art Reports (2014), pp. 21–
42. 1[JWJ∗ 14] J ONES B., WARD S., JALLEPALLI A., P ERENIA J.,
BARGTEIL A. W.: Deformation embedding for point-basedelastoplastic simulation. ACM Transactions on Graphics 33, 2(2014), 21:1–21:9. 3[LP12] L AIBE G., P RICE D. J.: Dusty gas with smoothed particle
hydrodynamics - II. implicit timestepping and astrophysical dragregimes. Monthly Notices of the Royal Astronomical Society 420,
3 (2012). 3
c 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.

[MCG03] M ÜLLER M., C HARYPAR D., G ROSS M.: Particlebased fluid simulation for interactive applications. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on
Computer Animation (2003), pp. 154–159. 2, 3, 5
[Mon89] M ONAGHAN J. J.: On the problem of penetration inparticle methods. Journal of Computational Physics 82, 1 (1989),
1–15. 2
[Mon92] M ONAGHAN J. J.: Smoothed particle hydrodynamics.
Annual Review of Astronomy and Astrophysics 30 (1992), 543–
574. 2, 3, 6
[Mon97] M ONAGHAN J. J.: Implicit SPH drag and dusty gas
dynamics. Journal of Computational Physics 138, 2 (Dec. 1997),
801–820. 3[MP89] M ILLER G., P EARCE A.: Globular dynamics: A connected particle system for animating viscous fluids. Computers& Graphics 13, 3 (1989), 305–309. 2[MSKG05] M ÜLLER M., S OLENTHALER B., K EISER R.,
G ROSS M.: Particle-based fluid-fluid interaction. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics Symposium on
Computer Animation (2005), pp. 237–244. 1, 2, 3
[OK12] O RTHMANN J., KOLB A.: Temporal blending for adaptive SPH. Computer Graphics Forum 31, 8 (2012), 2436–2449.
[PPLT06] PAIVA A., P ETRONETTO F., L EWINER T., TAVARES
G.: Particle-based non-newtonian fluid animation for melting
objects. In Proceedings of SIBGRAPI 2006. (2006), pp. 78–85.
2, 3[Raf07] An incompressible SPH method for simulation of unsteady viscoelastic free-surface flows. International Journal of
Non-Linear Mechanics 42, 10 (2007), 1210–1223. 3[REN∗ 04] R ASMUSSEN N., E NRIGHT D., N GUYEN D.,
M ARINO S., S UMNER N., G EIGER W., H OON S., F EDKIW R.:
Directable photorealistic liquids. In Proceedings of the 2004
ACM SIGGRAPH/Eurographics Symposium on Computer Animation (2004), pp. 193–202. 2, 3
[RLY∗ 14] R EN B., L I C., YAN X., L IN M. C., B ONET J., H U
S.-M.: Multiple-fluid SPH simulation using a mixture model.
ACM Transactions on Graphics 33, 5 (2014), 171:1–171:11.
[SG11] S OLENTHALER B., G ROSS M.: Two-scale particle simulation. ACM Transactions on Graphics 30, 4 (2011), 81:1–81:8.
1[SP08] S OLENTHALER B., PAJAROLA R.:
Density contrast SPH interfaces.
In Proceedings of the 2008 ACM
SIGGRAPH/Eurographics Symposium on Computer Animation(2008), pp. 211–218. 1[SP09] S OLENTHALER B., PAJAROLA R.: Predictive-corrective
incompressible SPH. ACM Transactions on Graphics 28, 3(2009), 40:1–40:6. 1, 3, 4
[SSJ∗ 14] S TOMAKHIN A., S CHROEDER C., J IANG C., C HAI L.,
T ERAN J., S ELLE A.: Augmented MPM for phase-change andvaried materials. ACM Transactions on Graphics 33, 4 (2014),
138:1–138:11. 2, 3, 4, 5
[SSP07] S OLENTHALER B., S CHLÄFLI J., PAJAROLA R.: A unified particle model for fluid-solid interactions. Computer Animation and Virtual Worlds 18, 1 (2007), 69–82. 1, 3
[TNF14] TAKAHASHI T., N ISHITA T., F UJISHIRO I.: Fast simulation of viscous fluids with elasticity and thermal conductivity using position-based dynamics. Computers & Graphics 43(2014), 21–30. 2[TPF91] T ERZOPOULOS D., P LATT J., F LEISCHER K.: Heating
and melting deformable models. Journal of Visualization andComputer Animation 2 (1991), 68–73. 2

Takahashi et al. / Implicit Formulation for SPH-based Viscous Fluids[WRK∗ 10] W ICKE M., R ITCHIE D., K LINGNER B. M., B URKE
S., S HEWCHUK J. R., O’B RIEN J. F.: Dynamic local remeshingfor elastoplastic simulation. ACM Transactions on Graphics 29,
4 (2010), 49:1–49:11. 2[WT08] W OJTAN C., T URK G.: Fast viscoelastic behavior withthin features. ACM Transactions on Graphics 27, 3 (2008), 47:1–
47:8. 2[ZLKW13] Z HOU Y., L UN Z., K ALOGERAKIS E., WANG R.:
Implicit integration for particle-based simulation of elasto-plastic
solids. Computer Graphics Forum 32, 7 (2013), 215–223. 3

Appendix A: Details of Coefficient Extraction
Due to complexity of extracting coefficients and similarity
of computing coefficients, we divide them into three groups:
Gi (cui ui , cvi ui , and cwi ui ), G j (cu j ui , cv j ui , and cw j ui ), and Gk(cuk ui , cvk ui , and cwk ui ).
By scanning particle j once from particle i, we can easily
compute and extract coefficients in Gi .
Since computing coefficients in G j requires ai j and ωi j(summation of ∇Wi j over particle j), and ∇Wi j and α jk(summation of Vk ∇W jk over particle k), we need to compute ωi j and αi j (we can access α jk from particle j if αi j
is computed in a previous loop) in advance using anotherloop unlike the case of Gi . Therefore, we first compute ωi jand αi j in the first loop and then use them to compute andextract coefficients in G j in the second loop. Computed ωi jand αi j can also be used to obtain coefficients in Gi .
To compute coefficients in Gk , we need to access particle k from particle i. However, since particle k is a neighbor of particle j, we cannot directly access particle k fromparticle i. Hence, we access particle k via particle j. In thatcase, however, we cannot take a sum of quantities computed between particles i and j at particle k (another scanfor particle j at particle k leads to quadruple loops, which
are costly and make coefficient extractions more complex).
Therefore, we decompose coefficients cuk ui , cvk ui , and cwk ui
without taking a sum over particle j and separately extract
coefficients at particle k. Assuming that we are accessing
particle k∗ via particle j∗ , we can write a part of coefficients
cuk ui | j∗ k∗ , cvk ui | j∗ k∗ , and cwk ui | j∗ k∗ for cuk ui , cvk ui , and cwk ui
from Eqs. (10), (11), and (12) as
cuk ui | j∗ k∗ = −m̂µ j∗(2∇Wi j∗ ,x a j∗ k∗ ,x + ∇Wi j∗ ,y a j∗ k∗ ,y + ∇Wi j∗ ,z a j∗ k∗ ,z ),cvk ui | j∗ k∗ = −m̂µ j∗ ∇Wi j∗ ,y a j∗ k∗ ,x ,cwk ui | j∗ k∗ = −m̂µ j∗ ∇Wi j∗ ,z a j∗ k∗ ,x .
By adding the part of coefficients above to the matrix at particle k while also scanning particle j from particle i, we canextract all coefficients of the matrix.
Adding coefficients to the matrix separately does work.
However, this approach is inefficient because we generally
use structures specialized for a sparse matrix, e.g., CSR and
Coordinate list, and adding values to such structures frequently is costly. To avoid this, we use auxiliary storage,
which is associated with particle i, to preserve coefficients

Particle i

Particle j

0

1

2

4

3

2

i ’s storage
Add parts ofcoefficients
0
0

0

3

1

2

1

4

2

3

3

5

3

4

4

1

Particle k

Figure 8: Illustration of computational flow for extracting
coefficients in Gk . Particle k is accessed from particle i via
particle j. At particle k, parts of coefficients in Gk are added
to a cell in particle i’s storage, whose id matches k’s id, as
shown by red arrows. Black arrows represent accessible particles and storage.
from particle k to i with k’s id. We add parts of coefficients
to i’s storage at particle k, grouping them based on k’s id
to minimize the number of adding coefficients to the matrix. Then, after scans over particles j and k from particle iare finished, we add coefficients in Gk to the matrix using
the storage. In our experiments, minimizing the number ofaccess to the matrix using the auxiliary storage accelerated
our coefficient extractions by a factor of 4.5 as compared toadding coefficients to the matrix directly at particle k. A flow
for coefficient extractions for Gk is illustrated in Figure 8.
We
also
perform
similarprocedures
explained
above
toextractcoefficients
cui vi , cvi vi , cwi vi , cu j vi , cv j vi , cw j vi , cuk vi , cvk vi , and cwk vi forvi and cui wi , cvi wi , cwi wi , cu j wi , cv j wi , cw j wi , cuk wi , cvk wi , andcwk wi for wi . We show our algorithm for coefficient
extraction in Algorithm 3.
Algorithm 3 Algorithm for coefficient extraction
1: initialize a matrix
2: for all fluid particle i do
3:compute µ̂i , ωi j and αi j
4: compute m̂
5: for all fluid particle i do
6:
initialize storage for uk , vk , and wk
7:
add cui ui , cvi ui , cwi ui , cui vi , cvi vi , cwi vi , cui wi , cvi wi , andcwi wi to the matrix
8:
for all fluid particle j do
9:compute ∇Wi j and ai j
10:
add cu j ui , cv j ui , cw j ui , cu j vi , cv j vi , cw j vi , cu j wi , cv j wi ,and cw j wi to the matrix
11:
for all fluid particle k do
12:compute a jk
13:
add cuk ui , cvk ui , cwk ui , cuk vi , cvk vi , cwk vi , cuk wi , cvk wi ,and cwk wi to the i’s storage with k’s id
14:
for all i’s storage do
15:
add cuk ui , cvk ui , cwk ui , cuk vi , cvk vi , cwk vi , cuk wi , cvk wi ,and cwk wi to the matrix using the storagec 2015 The Author(s)Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.