FleCSPH notes

Hyun Lim, Julien Loiseau

May 20, 2020

1 FleCSI
FleCSI is a compile-time configurable framework designed to support multi-physics application development. As
such, FleCSI provides a very general set of infrastructure design patterns that can be specialized and extended to
suit the needs of a broad variety of solver and data requirements. FleCSI currently supports multi-dimensional
mesh topology, geometry, and adjacency information, as well as n-dimensional hashed-tree data structures, graph
partitioning interfaces, and dependency closures.
FleCSI introduces a functional programming model with control, execution, and data abstractions that are con-
sistent both with MPI and with state-of-the-art, task-based runtimes such as Legion, HPX and Charm++. The
abstraction layer insulates developers from the underlying runtime, while allowing support for multiple runtime
systems including conventional models like asynchronous MPI.
The intent is to provide developers with a concrete set of user-friendly programming tools that can be used now,
while allowing flexibility in choosing runtime implementations and optimization that can be applied to future
architectures and runtimes.
FleCSI’s control and execution models provide formal nomenclature for describing poorly understood concepts such
as kernels and tasks. FleCSI’s data model provides a low-buy-in approach that makes it an attractive option for
many application projects, as developers are not locked into particular layouts or data structure representations.
FleCSI currently provides a parallel but not distributed implementation of binary, quad- and octree topologies.
Domain decomposition is implemented using space-filling curves, such as Morton ordering curve.
At current stage, FleCSI framework requires implementation of a driver and a specialization driver. The role of
the specialization driver is to provide the data model and its parallel distribution. Currently, FleCSI does not has
this feature fully implemented, so we provide it. The next step will be to incorporate it directly from FleCSPH
to FleCSI as we reach a good level of performance. The driver represents the general execution of the resolution
without worrying of the data locality and communications. As FleCSI is a code in development, its structure may
change in the future and we keep track of these changes in FleCSPH1 .

2 Domain and tree construction

In the previous section we described the FleCSI framework. This part gives details on the domain decomposition
using octree, tree construction and search algorithm.

2.1 Domain decomposition

In the current version of FleCSI the domain decomposition is done using the Morton ordering. It allows to describe,
sort and distribute particles based on a unique value – key.
In principle, various space-filling curves can be used, with their own advantages and flaws:
• Morton ordering: interlace the bits of X, Y and Z positions to create the key. It is very simple to compute,
but the curve suffers discontinuous jumps.
• Hilbert-Peano: interlace the bits, but also add rotations. In this ordering, data locality in memory guarantees
locality in space.
1 This note is authorized for unlimited release under LA-UR-17-25908

i
 • Other space-filling curves: hexagonal space filling curves, ...?
This first implementation is based on the Morton ordering which is used during several steps:
• The distribution part, to be able to split the particles between the processes providing a good locality in the
data.
• The tree construction and search.

2.2 Binary, Quad- and Octrees

3 General algorithm
The main distributed algorithm is presented in algorithm 1

Algorithm 1 Main algorithm

1: procedure specialization driver(input parameter file p)
2: Read parameter file p
3: Set simulation parameters
4: Distributed read input data files, specified in parameter file
5: Set additional parameters specified in input data files
6: while iterations do
7: Distribute the particles using distributed quick sort . Using Morton keys
8: Compute total range
9: Generate the local tree
10: Share branches
11: Compute ghosts particles
12: Update ghosts data
13: Do physics
14: Update ghosts data
15: Do physics
16: Periodic analysis and output
17: end while
18: end procedure

In the current version the 7 is based on a distributed quick sort algorithm. Each process sends a sample of its keys
to the master (or submaster for larger cases). We have set this to 256 Kb of key data per process, but it can be
increased for larger simulations. After receiving sample keys, the master determines the general ordering for all
the processes and shares the pivots. Then each process sorts its local keys and, in a global communication step,
the particles are distributed to the process on which they belong. The advantage is that it is a quick distribution
algorithm, but it can lead to bad load balancing.
• The ordering may not be perfect in terms of the number of particles per processes. But by changing the
amount of data exchanged to the master can lead to better affectation (?)
• The load balancing also depends on the number of neighbors of each particles. If a particle is located in a
poor area with large space between the particles this can lead to bad load balancing too.
After the sorting step the local tree can be created on each process. To be able to look for the ghosts and shared
particles we need to share some information with the neighbors on the tree. At line (10), the algorithm searches
for the neighboring branches, which are affecting the local particles. We compute the global bounding box of each
processes and based on this information each process can then compute the affecting branches to share from its
local tree. This new information is then added to the local tree by considering NON LOCAL particles. This data
structure does just contain the position and mass of the distant particle.
The branch sharing allows to compute the ghosts for this step. Each process performs a local search in the tree and
computes the required ghosts particles (the NON LOCAL bodies). Those data for shared and ghosts are stored

ii
and are use to share the complete particle information when 12 is invoked. As the ghosts data remain the same
within an iteration, the 12 can be used several times to update local information on remote particles.

4 I/O
Large-scale simulations require an efficient, parallel and distributed I/O. We base this first implementation on HDF5
file structure with H5Part and H5Hut. The I/O was developed at LANL and provides a simple interface to read
and write data in H5Part format. The first requirement is to allow the user to work directly with the Paraview
visualization tool and splash.

5 SPH Formulation in FleCSPH

Smoothed particle hydrodynamics (SPH) is an explicit numerical meshfree method that solves partial differential
equations (PDE) of hydrodynamics by discretizing the flow with a set of fluid elements called particle. The main
SPH formula to interpolate a quantity A(~r), which is specified by its values on a set of particles Ab ≡ A(~rb ), is as
follows (Rosswog 2009): X
A(~r) ' Vb Ab W (|~r − ~rb |, h) (1)
b∈Ω(~
r)

where W is a smoothing kernel, h is the smoothing length (hydro interaction range) at a position ~r, and Vb is a
volume element, usually Vb = mb /ρb .
In comparison with Eulerian methods, SPH has several advantages. It can easily adapt to complex geometries,
naturally handle low density regions, and does not require low-density floor to handle vacuum. It conserves mass by
construction, and can be easily made to conserve linear momentum, angular momentum, and energy up to roundoff.
Another advantage of using SPH is its exact advection of fluid properties. Furthermore, the same tree which is used
to find particle neighbors, can be employed for computing Newtonian gravitational forces.
Shortcomings of SPH are its convergence which is restricted to low-order, and high sensitivity to the initial particles
distribution. Also, SPH struggles with resolving turbulence-dominated flows and requires special care when handling
high gradients, such as shocks and stellar surface.
The starting point at a continuum limit is Euler ideal fluid equations in the Lagrangian formulation, expressing
conservation equations of mass and linear momentum:
dρ
= −ρ∇ · ~v , (2)
dt  
du P dρ
= , (3)
dt ρ2 dt
d~v ∇P
=− + ~g , (4)
dt ρ
(5)

where d/dt = ∂t + ~v · ∇ and ~g is a gravitational acceleration. The latter may be due to interparticle gravitational
interaction, an external gravitational field, or both.

iii
5.1 Basic formulation : Vanilla Ice SPH
In its simplest form, SPH discretization uses the volume element Vb = mb /ρb , a constant smoothing length h, and
artificial viscosity term Πab :
X
ρa = mb Wab , (6)
b
 
dua X Pa 1
= mb + Πab ~vab · ∇a Wab , (7)
dt ρ2a 2
b
 
d~va X Pa Pb
=− mb + + Π ab ∇a Wab + ~
ga , (8)
dt ρ2a ρ2b
b

where Wab = W (|~ra − ~rb |, h). The viscous stress tensor Πab may be defined in different ways; currently we have
adopted the following form:
−αc̄ab µab +βµ2ab
(
ρ̄ab for ~rab · ~vab < 0,
Πab = (9)
0 otherwise,
where the following quantities are defined:

h̄ab~rab · ~vab
µab = , (10)
|~rab |2 + h̄2ab
~rab = ~ra − ~rb , ~vab = ~va − ~vb , (11)
average speed of sound: c̄ab = (ca + cb )/2, (12)
average density: ρ̄ab = (ρa + ρb )/2, (13)
average smoothing length: h̄ab = (ha + hb )/2, (14)
(15)

In the basic formulation with a constant smoothing length, h̄ab ≡ h. The values of , α, and β control the strength
and application of artificial viscosity. Their default values:  = 0.01, α = 1.0, and β = 2.0.
The quantity c̄ab is a speed of sound, averaged between particles a and b: c̄ab = (ca + cb )/2. It is computed as usual:
s 
∂P
ca = , (16)
∂ρ S,a

where the partial derivative is taken under constraint of a constant entropy S, with thermodynamic conditions at
particle a. E.g., the following (Newton-Laplace) equation,
s
ΓPa
ca = (17)
ρa

holds both for polytropic and ideal fluid equations of state (see Section 7).
If ~g = 0, then the basic formulation, when integrated with a symplectic integrator, conserves energy, momentum
and angular momentum exactly (e.g. Rosswog 2009, Section 2.4).
Alternatively, FleCSPH features an implementation of the so-called thermokinetic formulation, in which the total
particle energy is evolved: ea = ua + 12 va2 . Corresponding discretized version of the energy equation reads,
 
dea X Pa~vb Pb~va ~va + ~vb ~
=− mb + + Π ab · ∇a Wab . (18)
dt ρ2a ρ2b 2
b

5.2 The SPH Equations with grad-h Term

The smoothing length h is not necessarily constant. Rather, h should be adapted in space and time whenever
densities and length scales vary by large amounts. There are many literatures about ways to adjust the smoothing

iv
length but here we adapt (See Sec.2.6 in Rosswog 2009 for more details)
 1/3
ma
ha = η (19)
ρa
where η should be chosen in the range between 1.2 and 1.5.
Note that the SPH equations in previous section were derived under assumption that the smoothing lengths are
constant. So, we need different set of equations. The SPH equations can be directly derive from variational principle
with fluid Lagrangian. Unlike vanilla SPH formulation, we take the changes of h into account. This introduces
additional term called grad-h term. This additional term increases the accuracy of SPH and the conservation
properties int the presence of varying h. We refer reader to see Sec.3.3 in Rosswog 2009 for whole detail derivations.
Below is summary of the SPH equations with grad-h term
X
ρa = mb Wab (rab , ha ) (20)
a

where smoothing length ha is

 1/3
ma
ha = η (21)
ρa
The energy equation is
dua 1 Pa X
= mb~vab · ∇a Wab (ha ) (22)
dt Ωa ρ2a
b

The momentum equation is

 
d~va X Pa Pb
=− mb ∇a Wab (ha ) + ∇a Wab (hb ) (23)
dt Ωa ρ2a Ωb ρ2b
b

where
∂ha X ∂Wab (ha )
Ωa ≡ 1 − mb (24)
∂ρa ∂ha
b

6 Kernels
There are many choices for selecting a smoothing kernel. In the code, the kernel choice is specified by a sph kernel
parameter, and the smoothing length is defined to be equal to the kernel support radius. Note that some works
define smoothing length differently, such that the kernel support radius becomes a multiple of smoothing length.
For example, for a cubic spline kernel it is twice the h. We do not make such distinction; all kernels below must
satisfy the following normalization condition:
ZZZ
W (~r, h)dD ~r = 1, (25)
SD (h)

where the integration is performed over the D-dimensional volume of a sphere SD (h) of radius h.
All the currently implemented kernels possess spherical symmetry, which makes it easy to impose exact conservation
of linear momentum. For such kernels, the gradients can be computed as follows:
∇a Wab ≡ ∇a W (|~ra − ~rb |, ha ) (26)
dW
= ~εab , (27)
dr
where ~εab ≡ ~rab /|~rab | is a unit vector in the direction from particle b to particle a. It is therefore sufficient to
implement dW/dr for each kernel.
In the formulae below, we define q ≡ |~r|/h.

v
6.1 sph kernel = "cubic spline"
The simplest (but not the best one) is a Monaghan’s cubic spline kernel:

1 − 6q 2 + 6q 3 if 0 ≤ q ≤ 1/2,
σD 
W (~r, h) = D 2(1 − q)3 if 1/2 ≤ q ≤ 1, (28)
h 
0 otherwise,


where q = r/h, D is the number of dimensions and σD is a normalization constant:

 
4 40 8
σD = , , in 1D, 2D and 3D resp. (29)
3 7π π

Radial derivative of the cubic spline kernel is:


dW σD −6q(2 − 3q) if 0 ≤ q ≤ 1/2,

= D+1 −6(1 − q)2 if 1/2 ≤ q ≤ 1, (30)
dr h 
0 otherwise.


6.2 sph kernel = "gaussian"

The Gaussian kernel is: ( 2
σD e−q if 0 ≤ q ≤ 3,
W (~r, h) = D (31)
h 0 if q > 3
where σD is :  
1 1 1
σD = 1/2
, , in 1D, 2D and 3D resp. (32)
π π π 3/2
Radial derivative of the cubic spline kernel is:
( 2
dW σD −2qe−q if 0 ≤ q ≤ 3,
= D+1 (33)
dr h 0 if q > 3

6.3 sph kernel = "quintic spline"

The Quintic spline kernel is:
 5 5 5
[(3 − q) − 6(2 − q) + 15(1 − q) ] if 0 ≤ q ≤ 1,

σD [(3 − q)5 − 6(2 − q)5 ]

if 1 ≤ q ≤ 2,
W (~r, h) = 5
(34)
hD (3 − q)
 if 2 ≤ q ≤ 3,

0 if q > 3


where σD is:  
1 7 3
σD = , , in 1D, 2D and 3D resp. (35)
120 478π 359π
Radial derivative of the cubic spline kernel is:


 [−5(3 − q)4 + 30(2 − q)4 − 75(1 − q)4 ] if 0 ≤ q ≤ 1,
[−5(3 − q)4 + 30(2 − q)4 ]

dW σD  if 1 ≤ q ≤ 2,
= D+1 4
(36)
dr h 
 −5(3 − q) if 2 ≤ q ≤ 3,

0 if q > 3


vi
6.4 sph kernel = "Wendland C2"
The Wendland C2-continuous kernel (C2) for 2D and 3D is
(
4
σD 1 − 2q (2q + 1) if 0 ≤ q ≤ 2,
W (~r, h) = D (37)
h 0 if q > 2
where σD is:  
7 21
σD = , in 2D and 3D resp. (38)
4π 16π
For 1D: (
q 3


σD 1− 2 (1.5q + 1) if 0 ≤ q ≤ 2,
W (~r, h) = (39)
h 0 if q > 2
where σD is:
5
σD = in 1D (40)
8

6.5 sph kernel = "Super Gaussian"

The super Gaussian Kernal is given by
( 2
e−q D


1 2 + 1 − q2 if 0 ≤ q ≤ 3,
W (~r, h) = . (41)
π D/2 hD 0 if q > 3

6.6 sph kernel = "Wendland C4"

The Wendland C4-continuous kernel (C4) for 2D and 3D is
(
6
σD 1 − 2q ( 35 2
12 q + 3q + 1) if 0 ≤ q ≤ 2,
W (~r, h) = D (42)
h 0 if q > 2
where σD is:  
9 495
σD = , in 2D and 3D resp. (43)
4π 256π
For 1D: (
q 5


σD 1− 2 (2q 2 + 2.5q + 1) if 0 ≤ q ≤ 2,
W (~r, h) = (44)
h 0 if q > 2
where σD is:
3
σD = in 1D (45)
4

6.7 sph kernel = "Wendland C6"

The Wendland C6-continuous kernel (C6) for 2D and 3D is
(
8
σD 1 − 2q (4q 3 + 6.25q 2 + 4q + 1) if 0 ≤ q ≤ 2,
W (~r, h) = D (46)
h 0 if q > 2
where σD is:  
78 1365
σD = , in 2D and 3D resp. (47)
28π 512π
For 1D: (
q 7
( 21 19 2

3
σD 1− 2 8 q + + 4 q + 3.5q + 1) if 0 ≤ q ≤ 2,
W (~r, h) = (48)
h 0 if q > 2
where σD is:
55
σD = in 1D (49)
64

vii
 Dimensions b0 b1 b2 b3
1D −1.5404568 × 10−2 3.6632876 × 10−1 −4.6519576 × 10−4 −7.3658324 × 10−2
2D 5.2245027 × 10−2 1.3090245 × 10−1 1.9358485 × 10−2 −6.1642906 × 10−3
3D 2.7012593 × 10−2 2.0510827 × 10−2 3.7451957 × 10−3 4.7013839 × 10−2

6.8 sph kernel = "Sinc"

The sinc Kernal is given by (
sincn π


σD,n 2 if 0 ≤ q ≤ 2,
W (~r, h) = D , (50)
h 0 if q > 2
where n is kernel index which is ranging from 3 to 12. and σD is

1/2 −1/2
b0 + b1 n + b2 n + b3 n
 for 1D,
−1
σD,n = b0 + b1 n + b2 n + b3 n −2 (51)
for 2D,

b0 + b1 n1/2 + b2 n + b3 n3/2 for 3D,


where the values of coefficients b0 , b1 , b2 , and b3 as a function of the dimensions are given by

7 Equations of State
To understand the inner property of stars, one needs to find the equation which describes the relation between the
pressure of matter and its density, temperature and other compositions such that

P = P (ρ, T, Ye , ...) (52)

First, we consider analytic equations of state that are relevant for binary neutron stars.

7.1 Ideal Gas

Ideal gas equation of state is
P (ρ, u) = (Γ − 1)ρu (53)
where Γ is the adiabatic index of the gas. For a monatomic gas, we set Γ = 5/3. For another test such as sod tube,
people use Γ = 1.4

7.2 Piecewise Polytrope

For more relevant simulation, we choose piecewise polytropic EOS.(Ideal gas EOS is still good for many simple
test cases like Sod tube) In our case, we assume constant entropy so that many thermodynamic situations can be
approximated as polytropes or piecewise functions made up of polytropes.
For neutron star case, we assume degenerated Fermi gas of neutrons then polytropic constant for a non-relativistic
degenerated neutron gas is
(3π 2 )2/3 ~2
K0 = 8/3
(54)
5mn
where mn is the mass of a proton and ~ is a Planck constant. For polytropic index, we set
5
γ0 = (55)
3
In the relativistic case,
5
γ1 = (56)
2

viii
Then, piecewise polytrope EOS is (
K0 ργ0 ifρ ≤ ρ0
P (ρ) = γ
K0 ρ0 0 γ1
(57)
γ
ρ0 1
ρ ifρ > ρ0

where ρ0 = 5 × 1014 g/cm3 . We can combine the piecewise polytropic EOS with ideal gas to attain an EOS valid at
both low and high densities. For more realistic studies, we need to consider different types of analytic EOSs such
as Maxwell-Boltzmann and Helmholtz EOSs. Also, we will put the functionality that can control tabulated EOS.

7.3 Zero Temperature Equations of State

Another interesting problem using SPH is the double white dwarf (DWD) simulations for studying possible pro-
genitors to type Ia supernovae. Here, we use zero temperature equations of state (ZTWD) as a variation of the self
consistent field technique. In ZTWD, the electron degeneracy pressure P varies with the mass density ρ according
to the relation h i
P = A x(2x2 − 3)(x2 + 1)1/2 + 3 sinh−1 x (58)

where the dimensionless parameter

 ρ 1/3
x≡ (59)
B
and the constant A and B are
πm4e c5
A≡ = 6.00288 × 1022 dynes cm−2 (60)
3h3
B 8πmp  me c 3
≡ = 9.81011 × 105 g cm−3 (61)
µe 3 h

8 Initial Data
Initial particle configurations are constructed using various methods, including regular cubic lattices, random par-
ticle distributions, or sequences of spherical shells. For a star, the density as a function of radius ρ(~r) can be found
by solving Lane-Emden equation, where we use polytropic index n = 1 (see Section 9). We pick particles with
variable mass, such that the mass of a particle is computed as:
ρ(~
ri ) 3N
mi = with nr = (62)
nr 4πR3
The smoothing length is set to a constant and uniform for all particles:
r
1 3NN
h= (63)
2 4πn
Here we choose NN , the average number of neighbors, to be 100.

9 Solving Lane-Emden Equation

Lane-Emden equation describes equilibrium configuration of a star. This is an equation which determines the star
density as a function of radius.
As we consider the star as a polytropic fluid, we use the equation of Lane-Emden which is a form of the Poisson
equation:

d2 θ 2 dθ
+ + θn = 0 (64)
dξ 2 ξ dξ
With ξ and θ two dimensionless variables. There is only exact solutions for a polytropic index n = 0.5, 1 and 2. In
our work we use a polytropic index of 1 which can correspond to a NS simulation.
For n = 1 the solution of equation 64 is:

ix
 sin(ξ)
θ(ξ) = (65)
ξ
We note ξ1 = π, the first value of ξ as θ(ξ) = 0. θ(ξ) is also defined as:
 ρ(ξ)  n1 ρ(ξ)
θ(ξ) = = (66)
ρc ρc
With ρc the internal density of the star and ρ the density at a determined radius. ξ is defined as:
s r
4πG (n−1)/n 2πG
ξ = Ar = ρc ×r = × r (for n = 1)
K(n + 1) K

With K a proportionality constant.

From the previous equations we can write the stellar radius R as:
r r
K(n + 1) (1−n)/2 K
R= ρc ξ1 = × ξ1 (67)
4πG 2πG
(We note that for n = 1 the radius does not depend of the central density.)
If, for example, we use dimensionless units as G = R = M = 1 (for the other results we use CGS with G =
6.674 × 10−8 cm3 g −1 s−2 ) We can compute K as:

R2 2πG
K= (68)
ξ12

N S1 N S2 N S3 N S4
Radius (cm) R=G=M =1 1500000 1400000 960000
K 0.636619 95598.00 83576.48 39156.94

Then we deduce the density function of r as :

r
sin(A × r) 2πG
ρ(ξ) = × ρc with A =
A×r K
As we know the total Mass M , the radius R and the gravitational constant G we can compute the central density
as:

M A3
ρc =
4π(sin(AR) − ARcos(AR))
Then we normalize the results to fit R = M = G = 1: K 0 = K/(R2 G), m0i = mi /M , h0i = hi /R, x~i 0 = x~i /R

9.1 Roche lobe problem

To create Hydrostatic Equilibrium Models we use a different equation of motion. This version use Roche Lobe:

d~
vi F~ Grav F~ Hydro v~i
= i + i + F~iRoche − (69)
dt mi mi trelax
With trelax ≤ tosc ∼ (Gρ)−1/2 and where F~iRoche is:

ˆ + µqyi ~yˆ − µzi ~zˆ

F~iRoche = µ(3 + q)xi ~x
0
With µ to be determined (for us µ = 0.069) and q = M M = 1 as the two polytropes have the same total mass. This
is apply to each star to get the equilibrium and the simulate the tidal effect.

x
9.2 Darwin problem
This is the way we use to generate the final simulation. The equation of motion for the relaxation is now:

d~
vi F~ Grav F~ Hydro v~i
= i + i + F~iRot − (70)
dt mi mi trelax

With trelax same as before and F~Rot defined by:

F~Rot = Ω2 (xi ~x̂ + yi ~ŷ) (71)

q
G(M +M 0 )
With Ω = a3 .
Or Lz = Qzz Ω and Qzz = i (x2i + yi2 ). At t = 0 we compute the total angular moment Lz which stay constant.
P
Using it during the relaxation we can compute Ω as: Ω = QLzzz
just by recomputing Qzz .
2
Here the scheme is in N but just for the relaxation step.
For this relaxation we use two stars generated as before, applying equation of motion 70. Using a as the distance
between the two polytropes (Here a = 2.9 for R = 1) and ~x̂ going for the center of the first to the second star, and
~ẑ is like the rotation vector.

10 Time Integration Scheme

Algorithm 2 Leapfrog time integration

1: read initial state: {r, v, u} ← {r 0 , v 0 , u0 }
2: compute rho(r), P(rho, u), cs(rho, P, u)
3: save v 1/2 : v12 = v
4: compute acceleration: a(P,rho,v12)
5: compute Dudt(P,rho,v,v12)
6: while iterations do
7: select timestep dt ← ∆t
8: v += a*dt/2 . velocity kick: v n+1/2 = v n + an ∆t
2

n ∆t
9: u += Dudt*dt/2 . internal energy kick: un+1/2 = un + du dt 2
10: update neighbors
11: r += v*dt . drift: rn+1 = rn + v n+1/2 ∆t
12: compute rho(r), P(rho, u), cs(rho, P, u) . uses un+1/2 !
13: save v n+1/2 : v12 = v
14: compute acceleration: a(P,rho,v12) . an+1 = a[P n+1 , ρn+1 , Πab (v n+1/2 )]
15: v += a*dt/2 . second velocity kick: v n+1 = v n+1/2 + an+1 ∆t 2
16: update neighbors
du n+1


17: compute Dudt(P,rho,v,v12) . = u̇[P n+1 , ρn+1 , v n+1 , Πab (v n+1/2 )]
dt

n+1 ∆t
18: u += Dudt*dt/2 . second internal energy kick: un+1 = un+1/2 + du dt 2
19: end while

For the time integrator, we are using the leapfrog algorithm, specifically its ”kick-drift-kick” variation (see Algo-
rithm 2). This algorithm belongs to the family of symplectic integrators and, in the absence of gravity, conserves
energy exactly.
Leapfrog time integrators are efficient for Hamiltonian systems, specifically for particle simulations. The name
comes from the fact that particle velocities are updated at half-steps while the positions at integer steps, so that
the two leap over each other. A pair of updates from timestep n to n + 1 has the following form:

v n+1/2 = v n−1/2 + a(rn )∆t, (72)

n+1 n n+1/2
r =r +v ∆t. (73)

xi
Note that in this simple form the accelerations a(rn ) are computed synchronously with positions and assumed
independent from the velocities v n . This formulation is time-symmetric and reversible up to roundoff.
An equivalent ”kick-drift-kick” formulation was shown to be stable for variable time steps:

v n+1/2 = v n + a(rn )∆t/2, ”kick” (74)

n+1 n n+1/2
r =r +v ∆t, ”drift” (75)
n+1 n n+1
v = v + a(r )∆t/2, ”kick” (76)

In application to the basic SPH formulation (7-8), this method reads:

v n+1/2 = v n + a[rn , Πab (v n−1/2 )]∆t/2, (77)

du h n n i
un+1/2 = un + r , v , Πab (v n−1/2 ) ∆t/2, (78)
dt
rn+1 = rn + v n+1/2 ∆t, (79)
h i
v n+1 = v n+1/2 + a rn+1 , Πab (v n+1/2 ) ∆t/2, (80)
du h n+1 n+1 i
un+1 = un+1/2 + r ,v , Πab (v n+1/2 ) ∆t/2. (81)
dt
where we highlighted variable dependencies between different time levels.
In this formulation, the internal energy is stored and updated synchronously with the velocities rather than the
coordinates. This is done to achieve exact energy conservation:
X h i X1  2  2 
n+1/2 n−1/2 n+1/2 n−1/2 n+1/2 n−1/2
E −E = ma ua − ua + ma ~va − ~va (82)
a a
2
 n n−1/2 n−1/2
d~v n

X dua ~va + ~va
= ma + ~van · a , where ~van := . (83)
a
dt dt 2

If we substitute corresponding time derivatives from (7-8), this expression vanishes (as shown in e.g. Rosswog 2009,
their Section 2.4):
~van + ~vbn
 
X P a n Pb n
E n+1/2 − E n−1/2 = − ma mb ~
v + ~
v + Π ab
n
· ∇a Wab = 0. (84)
ρ2a b ρ2b a 2
a,b

For this expression to vanish, the viscosity tensors Πa b in (77) and (78) above (or in (80) and (80)) should be
identical. Even though they might depend on the values of internal energy at previous half-step, un−1/2 , it does not
affect tensor symmetry and energy conservation at half-steps. Similarly, the pressure in (84) depends on the density
at the current integer step and internal energy at previous half-step, but it does not violate energy conservation
as long as the same value of the pressure is used to compute accelerations and time derivatives of internal energy.
n n
At the same time, when computing derivaties of internal energy using expression (7), the dot products ~vab · ∇a Wab
n
must be computed with velocities ~va at integer timesteps, as in (83).

10.1 Adaptive timestep

The timestep can be adaptive and determined by ∆t = λCFL min(∆t1 , ∆t2 ), where λCFL ≈ 0.1 is a Courant factor
limit, and the timescales ∆t1 , ∆t2 are:
 
ha
∆t1 = min , (85)
a ca (1 + 1.2α) + 1.2 β maxb µab
s
ha
∆t2 = min , (86)
a ˙~v
| a|

where ca is a sound speed, α and β are viscosity parameters, and µab is a viscosity function, as defined in equation
(10).

xii
11 Computing Gravitational Forces
Gravitational force acting on a particle F~aGrav due to attraction by other particles is described by the Newton’s
formula, with G being the Newton’s gravitational constant:
X ma mb
F~aGrav = G ~rab . (87)
(|~ra − ~rb |)3
b

To avoid O(N 2 ) computational complexity, we use the Fast Multipole Method (FMM), as described below.

11.1 Fast Multipole Method

See Algorithm 3 for the basic outline of the method. In this algorithm, the macangle is the angle of the Multipole
Acceptance Criterion. The acceleration is, for a center of mass c, the sum of contributions from the local particles
and the distant cells with:
X mp .(r~c − r~p )
f~c (r~c ) = − (88)
p
|r~c − r~p |3

With p the particle inside this cell and cell the cells that are accepted with the MAC. Here we directly consider the
gravitational acceleration, we don’t take in account the mass of the center of mass c. And G = 1 in our context.
The acceleration at a point from this center of mass is based on taylor series:

∂ f~c 1 ∂ f~c
f~(~r) = f~c (r~c ) + || || · (~r − r~c ) + (~r − r~c )| · || || · (~r − r~c ) (89)
∂ r~c 2 ∂ r~c ∂ r~c
~
The Jacobi matrix || ∂∂~
fc
r || is then:
 
3(xc −xp )(xc −xp ) 3(yc −yp )(xc −xp ) 3(zc −zp )(xc −xp )
1− |rc −rp |2 − |rc −rp |2 − |r~c −r~p |2
X mp  − 3(xc −xp )(yc2−yp )
 3(y −y )(y −y ) 3(z −z )(y −y ) 
− |r~c −r~p | 1 − c |r~c p−r~p |c2 p − c |r~cp−r~p |c2 p  (90)
p
|r~c − r~p |3  3(xc −x p )(zc −zp ) 3(y −yp )(zc −zp ) 3(z −z )(z −z )

− |r~c −r~p |2 − c |r~c − r~p |2 1 − c |r~c p−r~p |c2 p

∂fci X mp h 3.(rci − rpi )(rcj − rpj ) i

|| j
|| = − δij − (91)
∂rc p
|r~c − r~p |3 |r~c − r~p |2

With δij the identity matrix with δij = 1 if i = j where i, j runs spatial index from 1 to 3. For example, r1 = x,
r2 = y, and r3 = z as usual sense. (We do not consider covariant form of this because we are not considering
spacetime).
~
The Hessian matrix || ∂ r~∂cf∂cr~c || is then:

5(xc −xp )3 5(xc −xp )2 (yc −yp ) 5(xc −xp )2 (zc −zp )
 
|r~c −r~p |2 − 3(xc − xp ) |r~c −r~p |2 − 3(yc − yp ) |r~c −r~p |2 − 3(zc − zp )
∂ 2 fcx X 3mp  5(x −x 2
p ) (yc −yp ) 5(xc −xp )(yc −yp )2 5(xc −xp )(yc −yp )(zc −zp )

|| || = −  c
| − r~p |2 − 3(yc − yp ) |r~c −r~p |2 − 3(xc − xp ) |r~c −r~p |2

∂rci ∂rcj p
|r~c − r~p | 5  r
~ c
5(xc −xp )2 (zc −zp ) 5(xc −xp )(yc −yp )(zc −zp ) 5(xc −xp )(zc −zp )2

|r~c −r~p |2 − 3(zc − zp ) |r~c −r~p |2 |r~c −r~p |2 − 3(x c − x p )
(92)

2
c −xp ) (yc −yp ) 5(xc −xp )(yc −yp )2 5(xc −xp )(yc −yp )(zc −zp )
 5(x 
|r~c −r~p |2 − 3(yc − yp ) |r~c −r~p |2 − 3(xc − xp ) |r~c −r~p |2
∂ 2 fcy X 3mp  5(xc −xp )(yc −zp )2 5(yc −yp )3 5(yc −yp )2 (zc −zp )

|| || = − 
|r~c −r~p |2 − 3(xc − xp ) |r~c −r~p |2 − 3(yc − yp ) |r~c −r~p |2 − 3(zc − zp )
∂rci ∂rcj p
|r~c − r~p |5 
5(xc −xp )(yc −yp )(zc −zp ) 5(yc −yp )2 (zc −zp ) 5(yc −yp )(zc −zp )2

|r~c −r~p |2 |r~c −r~p |2 − 3(zc − zp ) |r~c −r~p |2 − 3(yc − yp )
(93)

xiii
 2
c −xp ) (zc −zp ) 5(xc −xp )(yc −yp )(zc −zp ) 5(xc −xp )(zc −zp )2
 5(x 
|r~c −r~p |2 − 3(zc − zp ) |r~c −r~p |2 |r~c −r~p |2 − 3(xc − xp )
∂fcz X 3mp  5(xc −xp )(yc −zp )(zc −zp ) 5(yc −yp )2 (zc −zp ) 5(yc −yp )(zc −zp )2

|| j
|| = − 
|r~c −r~p |2 |r~c −r~p |2 − 3(zr − zp ) |r~c −r~p |2 − 3(yc − yp ) 
∂rci ∂rc |r~c − r~p |5  
p 5(xc −xp )(zc −zp )2 5(yc −yp )(zc −zp )2 5(zc −zp )3
|r~c −r~p |2 − 3(xc − xp ) |r~c −r~p |2 − 3(yc − yp ) |r~c −r~p | 2 − 3(zc − zp )
(94)

" #
∂ 2 fci X 3mp 5(rci − rpi )(rcj − rpj )(rck − rpk ) 3
δij (rck − rpk ) + δjk (rci − rpi ) + δik (rcj − rpj )


|| || = − − (95)
∂rcj ∂rck p
|r~c − r~p |5 |r~c − r~p |2 w

where w = δij + δjk + δjk + ijk and ijk is 3D Levi-Civita symbol. We add Levi-Civita symbol to avoid zero in
denominator Here again, latin indices i, j, and k indicates spatial components

12 Resolution order
The resolution is done in this order:

Load data from file

while TotalTime not reached do
Move particles
Apply rotation with defined angular velocity
Compute density ρ
Compute pressure P and sound speed c
Compute F~hydro
Compute F~grav
if relaxation then
Compute F~roche or F~rot
end if
Compute acceleration ~a with the equation of motion
Compute ∆t
Compute new velocity from new acceleration
if output step then
output data to file
end if
TotalTime ← TotalTime + ∆t
end while

13 Applications
13.1 Sod Shock Tube
The Sod shock tube is the test consists of a one-dimensional Riemann problem with the following initial parameters
(
(1.0, 0.0, 1.0) if0 < x ≤ 0.5
(ρ, v, p)t=0 = (96)
(0.125, 0.0, 0.1) if0.5 < x < 1.0

This link shows some references and values for Sod shock tube problem that also includes boundary and jump
conditions([Link]
Also, we would like to re-generate the shock test result from Rosswog’s paper. In that paper, he shows the result of
a 2D relativistic shock tube test where the left state is given by [P, vx , vy , N ]L = [40/3, 0, 0, 10] and the right state
by [P, vx , vy , N ]R = [10−6 , 0, 0, 1] with Γ = 5/3
In our code, we use below parameters to get results

xiv
Algorithm 3 Gravitation computation
1: procedure tree traversal grav(branch sink)
2: if [Link] < Mcellmax then . Another choice criterion can be use
3: f~c ← ~0
∂ f~c ~
4: ∂~r ←0
~ ∂ f~c
5: TREE TRAVERSAL C2C(sink,[Link],f~c , ∂ fc ) ∂~
r . Compute f~c and ∂~ r using MAC
f~c
6: SINK TRAVERSAL C2P(sink, f~c , ∂∂~r ) . Expand to the particles below
7: else
8: for All children c of sink do
9: TREE TRAVERSAL GRAV(c)
10: end for
11: end if
12: end procedure
13:
14: function MAC(branch sink,branch source,double macangle)
15: dmax ← [Link] × 2
16: dist ← distance([Link], [Link])
17: return dmax /dist < macangle
18: end function
19:
20: procedure tree traversal c2c(branch sink, branch source, acceleration f~c )
21: if MAC(sink, source, macangle) then
22: f~c ← f~c + (− [Link]×([Link]−[Link])
|[Link]−[Link]|3 )
∂ f~c
23: ∂~ r← ...
24: else
25: if [Link] leaf () then
26: for All particles p of source do
27: f~c ← f~c + (− [Link]×([Link]−[Link])
|[Link]−[Link]|3 )
∂ f~c
28: ∂~ ← ...
r
29: end for
30: else
31: for All children c of source do
32: TREE TRAVERSAL C2C(sink, c, f~c )
33: end for
34: end if
35: end if
36: end procedure
37:
38: procedure TREE TRAVERSAL C2P(branch current, acceleration f~c )
39: if [Link] leaf () then
40: for All particle p of current do
δ f~c
41: [Link] ← f~c + δ[Link] .([Link] − [Link]) + ...
42: end for
43: else
44: for All children c of current do
~
45: TREE TRAVERSAL C2P(c, f~c , ∂∂ fr~cc )
46: end for
47: end if
48: end procedure

xv
13.2 Sedov Blast Wave
A blast wave is the pressure and flow resulting from the deposition of a large amount of energy in a small very
localized volume. This is another great test problem for computational fluid dynamics field.
There are different version of blast wave test but we consider the analytic solution for a point explosion is given by
Sedov, making the assumption that the atmospheric pressure relative to the pressure insider the explosion negligible.
The position of the shock as a function of time t, relative to the initiation of the explosion, is given by
1
et2
  δ+2
R(t) = (97)
ρ0

with δ = 2 and δ = 3 for cylindrical and spherical geometry respectively. The initial density ρ0 whereas e is a
dimensionless energy. Right behind the shock we ahve the following properties
Γ+1 2 2
ρ2 = ρ0 P 2 = ρ 0 w 2 v2 = w (98)
Γ−1 Γ+1 Γ+1
where the shock velocity is
dR 2 R(t)
w(t) = = (99)
dt δ+2 t
In numerical simulations, energy deposition in a single point is difficult to achieve. A solution to the problem is to
make use of the bursting balloon analogue. Rather than depositing the total energy in a single point, the energy is
released into a balloon of finite volume V
(P − P0 )V
e= (100)
Γ−1
The energy release in a balloon of radius r0 raises the pressure to the value

3(Γ − 1)e
P = (101)
(δ + 1)πr0δ

Here, we test 2D blast wave test. In this simulation, we use ideal gas EOS with Γ = 5/3 and we are assuming that
the undistributed area is at rest with a pressure P0 = 1.0−5 . The density is constant ρ0 , also in the pressurized
region.

xvi