Computational Astrophysics - Solving for Gravity

Computational Astrophysics - Solving for Gravity

Alexander Knebe, Universidad Autonoma de Madrid

1. Computational Astrophysics
Solving for GravityAlexander Knebe, Universidad Autonoma de Madrid

2. Computational Astrophysics
Solving for Gravity
•collisionless matter (e.g. dark matter)
•collisional matter (e.g. gas)
€
∂ρE()
∂t

€
+∇⋅ρE+p+
1
2µ
B
2
%
&
'
(
)
*
!
v −
1
µ
!
v ⋅
!
B
[]
!
B
,
-
.
/
0
1

€
∂ρ
!
v ()
∂t
€
∂ρ
∂t
€
+∇⋅ρ
!
v ⊗
!
v +p+
1
2µ
B
2
&
'
(
)
*
+
"
1 −
1
µ
!
B ⊗
!
B
&
'
(
)
*
+
€
+∇⋅ρ
!
v ()
€
=ρ−∇φ()
€
=0

€
=ρ
!
v ⋅−∇φ()+Γ−L()

€
d
!
x
DM
dt
=
!
v
DM

€
d
!
v
DM
dt
=−∇φ
!full set of equations
the equations
€
p=γ−1()ρε
€
ρε=ρE−
1
2
ρv
2
€
=4πGρ
tot
€
Δφ

€
∂
!
B
∂t
=−∇×
!
v ×
!
B
()
•Poisson’s equation
•Maxwell’s equation
•ideal gas equations

3. Computational Astrophysics
Solving for Gravity
•collisionless matter (e.g. dark matter)
•collisional matter (e.g. gas)
€
∂ρE()
∂t

€
+∇⋅ρE+p+
1
2µ
B
2
%
&
'
(
)
*
!
v −
1
µ
!
v ⋅
!
B
[]
!
B
,
-
.
/
0
1

€
∂ρ
!
v ()
∂t
€
∂ρ
∂t
€
+∇⋅ρ
!
v ⊗
!
v +p+
1
2µ
B
2
&
'
(
)
*
+
"
1 −
1
µ
!
B ⊗
!
B
&
'
(
)
*
+
€
+∇⋅ρ
!
v ()
€
=ρ−∇φ()
€
=0

€
=ρ
!
v ⋅−∇φ()+Γ−L()

€
d
!
x
DM
dt
=
!
v
DM

€
d
!
v
DM
dt
=−∇φ
!full set of equations
the equations
€
p=γ−1()ρε
€
ρε=ρE−
1
2
ρv
2
€
=4πGρ
tot
€
Δφ

€
∂
!
B
∂t
=−∇×
!
v ×
!
B
()
•Poisson’s equation
•Maxwell’s equation
•ideal gas equations
…and the force

4. Computational Astrophysics
Solving for Gravity

€
!
F (
!
x )=−m∇Φ(
!
x )
ΔΦ(
!
x )=4πGρ(
!
x )
Poisson’s Equation
!Poisson’s equation

5. Computational Astrophysics
Solving for Gravity

€
!
F (
!
x )=−m∇Φ(
!
x )
ΔΦ(
!
x )=4πGρ(
!
x )
Poisson’s Equation
grid approach

€
ΔΦ(
!
x
i,j,k
)=4πGρ(
!
x
i,j,k
)
!
F (
!
x
i,j,k
)=−m∇Φ(
!
x
i,j,k
)

€
(
!
x
i,j,k
=position of centre of grid cell (i,j,k))

€
!
F (
!
x
i
)=−
Gm
i
m
j
(x
i
−x
j
)
3
(
!
x
i
−
!
x
j
)
i≠j
∑
particle approach
!Poisson’s equation

6. Computational Astrophysics
Solving for Gravity

€
!
F (
!
x )=−m∇Φ(
!
x )
ΔΦ(
!
x )=4πGρ(
!
x )
Poisson’s Equation
grid approach

€
ΔΦ(
!
x
i,j,k
)=4πGρ(
!
x
i,j,k
)
!
F (
!
x
i,j,k
)=−m∇Φ(
!
x
i,j,k
)

€
(
!
x
i,j,k
=position of centre of grid cell (i,j,k))

€
!
F (
!
x
i
)=−
Gm
i
m
j
(x
i
−x
j
)
3
(
!
x
i
−
!
x
j
)
i≠j
∑
particle approach
!Poisson’s equation
weapon of choice: tree codes

7. Computational Astrophysics
Solving for Gravity

€
!
F (
!
x )=−m∇Φ(
!
x )
ΔΦ(
!
x )=4πGρ(
!
x )
Poisson’s Equation
grid approach

€
ΔΦ(
!
x
i,j,k
)=4πGρ(
!
x
i,j,k
)
!
F (
!
x
i,j,k
)=−m∇Φ(
!
x
i,j,k
)

€
(
!
x
i,j,k
=position of centre of grid cell (i,j,k))

€
!
F (
!
x
i
)=−
Gm
i
m
j
(x
i
−x
j
)
3
(
!
x
i
−
!
x
j
)
i≠j
∑
particle approach
!Poisson’s equation
weapon of choice: AMR codes

8. Computational Astrophysics
Solving for Gravity
the particle-mesh (PM) method

9. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)

10. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)

11. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)∀i∈N

12. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)
€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid
∀i∈N

13. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)
€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid

€
!
F (
!
g
k,l,m
)
3. differentiate potential to get forces
∀i∈N

14. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)
€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid

€
!
F (
!
g
k,l,m
)
3. differentiate potential to get forces

€
!
F (
!
g
k,l,m
)→
!
F (
!
x
i
)
4. interpolate forces back to particles
∀i∈N
∀i∈N

15. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation

€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid

€
!
F (
!
g
k,l,m
)
3. differentiate potential to get forces

€
!
F (
!
g
k,l,m
)→
!
F (
!
x
i
)
4. interpolate forces back to particles
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)
sounds like a waste of time and computer resources,
but exceptionally fastin practice
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)

16. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation

€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid

€
!
F (
!
g
k,l,m
)
3. differentiate potential to get forces

€
!
F (
!
g
k,l,m
)→
!
F (
!
x
i
)
4. interpolate forces back to particles
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)

17. Computational Astrophysics
Solving for Gravity
xaxis
!density assignment schemes

€
!
x
i
→ρ(
!
g
k,l,m
)
m
€
M(g
k
)=mW(d)d=x−g
k
ρ(g
k
)=
M(g
k
)
H
H
example:1 particle on 1 dimensional grid
Particle-Mesh Method
!
particle position xgrid point gk

18. Computational Astrophysics
Solving for Gravity
xaxis
!density assignment schemes

€
!
x
i
→ρ(
!
g
k,l,m
)
m
€
M(g
k
)=mW(d)d=x−g
k
ρ(g
k
)=
M(g
k
)
H
H
mass assignment function
example:1 particle on 1 dimensional grid
Particle-Mesh Method
!
particle position xgrid point gk

19. Computational Astrophysics
Solving for Gravity
xaxis
!density assignment schemes

€
!
x
i
→ρ(
!
g
k,l,m
)
m
H
example:1 particle on 1 dimensional grid
Particle-Mesh Method
•hierarchy of mass assignment schemes:
!Nearest-Grid-PointNGP
!Could-In-CellCIC
!Triangular-Shaped CloudTSC
!…
particle position xgrid point gk
"

20. Computational Astrophysics
Solving for Gravity
xaxis
!density assignment schemes
Nearest-Grid-Point (NGP):
m
H

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method
mass assignment function:

€
W(d)=
1d≤H/2
0otherwise
#
$
%

21. Computational Astrophysics
Solving for Gravity
xaxis
!density assignment schemes
Nearest-Grid-Point (NGP):
m
H

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method
particle shape:
€
S(x)=δ(x)
mass assignment function:

€
W(d)=
1d≤H/2
0otherwise
#
$
%

22. Computational Astrophysics
Solving for Gravity
xaxis
!density assignment schemes
Cloud-In-Cell (CIC):
m
H

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method
mass assignment function:

€
W(d)=
1−
d
H
d≤H
0 otherwise
$
%
&
&
'
&
&

23. Computational Astrophysics
Solving for Gravity
particle shape:

€
S(x)=
1x≤H/2
0otherwise
#
$
%
xaxis
!density assignment schemes
Cloud-In-Cell (CIC):
m
H

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method
mass assignment function:

€
W(d)=
1−
d
H
d≤H
0 otherwise
$
%
&
&
'
&
&

24. Computational Astrophysics
Solving for Gravity
xaxis
!density assignment schemes
Triangular-Shaped-Cloud (TSC):
m
H

€
W(d)=
3
4
−
d
H
#
$
%
&
'
(
2
d≤
H
2
1
2
3
2
−
d
H
#
$
%
&
'
(
2
H
2
≤d≤
3H
2
0 otherwise
*
+
,
,
,
,
-
,
,
,
,

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method
mass assignment function:

25. Computational Astrophysics
Solving for Gravity
particle shape:

€
S(x)=
1−
x
H
x≤H
0otherwise
$
%
&
'
&
xaxis
!density assignment schemes
Triangular-Shaped-Cloud (TSC):
m
H

€
W(d)=
3
4
−
d
H
#
$
%
&
'
(
2
d≤
H
2
1
2
3
2
−
d
H
#
$
%
&
'
(
2
H
2
≤d≤
3H
2
0 otherwise
*
+
,
,
,
,
-
,
,
,
,

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method
mass assignment function:

26. Computational Astrophysics
Solving for Gravity
!density assignment schemes

€
!
d =
!
x
i
−
!
g
k,l,m
M(
!
g
k,l,m
)=m
i
W(|d
x
|)W(|d
y
|)W(|d
z
|)
i=1
N
∑
ρ(
!
g
k,l,m
)=
M(
!
g
k,l,m
)
H
3
Nparticles on 3 dimensional grid

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method

27. Computational Astrophysics
Solving for Gravity
!density assignment schemes

€
!
d =
!
x
i
−
!
g
k,l,m
M(
!
g
k,l,m
)=m
i
W(|d
x
|)W(|d
y
|)W(|d
z
|)
i=1
N
∑
ρ(
!
g
k,l,m
)=
M(
!
g
k,l,m
)
H
3
Nparticles on 3 dimensional grid

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method
for every grid point we need to loop over all Nparticles…

28. Computational Astrophysics
Solving for Gravity
!density assignment schemes -in practice
Particle-Mesh Method

€
!
x
i
→ρ(
!
g
k,l,m
)
…rather loop over all particles
and
assign them to the appropriate grid points,
because the mapping xi →gkis rather easy

29. Computational Astrophysics
Solving for Gravity

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method
m
!!
!!
!density assignment schemes -in practice
•example for CIC assignment in 2D:

€
!
x
i
contributes its mass m
i
to the 4 closest grid points:

30. Computational Astrophysics
Solving for Gravity
!density assignment schemes -in practice
•example for CIC assignment in 3D:
Particle-Mesh Method

€
!
x
i
→ρ(
!
g
k,l,m
)

€
!
x
i
contributes its mass m
i
to the 8 closest grid points:

€
!
g
k,l,m+1
!
g
k+1,l,m+1
!
g
k,l,m
!
g
k+1,l,m

€
!
g
k,l+1,m+1
!
g
k+1,l+1,m+1
!
g
k,l+1,m
!
g
k+1,l+1,m

31. Computational Astrophysics
Solving for Gravity
!density assignment schemes
•which scheme to choose?
NGP = stepwise force(1 grid point)
CIC = continuous piecewise linear force(8 grid points)
TSC = continuous force and first derivative(27 grid points)
Particle-Mesh Method

32. Computational Astrophysics
Solving for Gravity
!density assignment schemes
•which scheme to choose?
NGP = too crude
CIC = common choice
TSC = pretty smooth
increased smoothing of density field
Particle-Mesh Method

33. Computational Astrophysics
Solving for Gravity
!density assignment schemes
•which scheme to choose?
NGP = too crude
CIC = common choice
TSC = pretty smooth
smoothing the density field will lead to a “bias” in the forces
but at the same time decrease the “variance”!
increased smoothing of density field
Particle-Mesh Method
bias=!F(!x)−!Ftrue(!x)( )∝εαvar=!F2(!x)−!F(!x)2∝N−β

34. Computational Astrophysics
Solving for Gravity
!density assignment schemes
•which scheme to choose?
NGP = too crude
CIC = common choice
TSC = pretty smooth
Particle-Mesh Method
(interplay between N and ! : N!3=const.)
smoothing the density field will lead to a “bias” in the forces
but at the same time decrease the “variance”!
bias=!F(!x)−!Ftrue(!x)( )∝εαvar=!F2(!x)−!F(!x)2∝N−β

35. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation

€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid

€
!
F (
!
g
k,l,m
)
3. differentiate potential to get forces

€
!
F (
!
g
k,l,m
)→
!
F (
!
x
i
)
4. interpolate forces back to particles
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)

36. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
€
ΔΦ
k,l,m
=ρ
k,l,m
Particle-Mesh Method
•relaxation technique: applicable and usable for anydifferential equation
•FTT technique: only applicable and usable for linear differential equation

37. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
€
ΔΦ
k,l,m
=ρ
k,l,m
Particle-Mesh Method
•relaxation technique: applicable and usable for anydifferential equation
•FTT technique: only applicable and usable for linear differential equation

38. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
€
ΔΦ
k,l,m
=ρ
k,l,m
Particle-Mesh Method
•Green’s function method:
-solve differential equation by Fourier transformation
-applicable and usable for linear differential equations

39. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation fast fourier transform method
•Green’s function method
Particle-Mesh Method
€
ΔΦ=ρ

€
Φ(
!
x )=G(
!
x −
!
$ x )∫∫∫ ρ(
!
$ x )d
3
$ x
;G(!x)=14πx

40. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation fast fourier transform method
•Green’s function method
€
ΔΦ=ρ
Particle-Mesh Method

€
Φ(
!
x )=G(
!
x −
!
$ x )∫∫∫ ρ(
!
$ x )d
3
$ x
;G(!x)=14πx

€
Φ=ρ⊗G
FFT "convolution becomes multiplication

€
ˆ
Φ =ˆ ρ
ˆ
G

€
ˆ
G =−
1
k
2
(with for Poisson’s equation )

41. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation fast fourier transform method
•Green’s function method
€
ΔΦ=ρ
Particle-Mesh Method

€
Φ(
!
x )=G(
!
x −
!
$ x )∫∫∫ ρ(
!
$ x )d
3
$ x
;G(!x)=14πx

€
Φ=ρ⊗G
FFT "convolution becomes multiplication

€
ˆ
Φ =ˆ ρ
ˆ
G

€
ˆ
G =−
1
k
2
(with for Poisson’s equation )
FFT-1
€
Φ
(FFT demands a regulargrid though!)

42. Computational Astrophysics
Solving for Gravity

€
ˆ
G (
!
g
k,l,m
)=−
1
sin
2k
x
2
#
$
%
&
'
( +sin
2
k
y
2
#
$
%
&
'
( +sin
2k
z
2
#
$
%
&
'
(

€
ˆ
G (
!
g
0,0,0
)=0,k
x
=
2πk
L
,k
y
=
2πl
L
,k
y
=
2πm
L

€
ˆ
G (
!
k )=−
1
k
2
€
" → "
Particle-Mesh Method
!numerically integrate Poisson’s equation fast fouriertransform method
•discretized Green’s function
compensates for grid anisotropies...

43. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
€
ΔΦ
k,l,m
=ρ
k,l,m
Particle-Mesh Method
•relaxation technique: applicable and usable for anydifferential equation
•FTT technique: only applicable and usable for linear differential equation

44. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
€
ΔΦ
k,l,m
=ρ
k,l,m
Particle-Mesh Method
obtain iterative solver by discretizing differential equation
!numerically integrate Poisson’s equation relaxation technique

45. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation relaxation technique
€
ΔΦ
k,l,m
=ρ
k,l,m
Particle-Mesh Method
€
ΔΦ
k,l,m
=∇⋅∇Φ
k,l,m
=
∂
∂x
∂
∂y
∂
∂z
'
(
)
)
)
)
)
)
*
+
,
,
,
,
,
,
⋅
∂Φ
k,l,m
∂x
∂Φ
k,l,m
∂y
∂Φ
k,l,m
∂z
'
(
)
)
)
)
)
)
*
+
,
,
,
,
,
,
=
1
H
∂
∂x
∂
∂y
∂
∂z
'
(
)
)
)
)
)
)
*
+
,
,
,
,
,
,
⋅
Φ
k+
1
2
,l,m
−Φ
k−
1
2
,l,m
Φ
k,l+
1
2
,m
−Φ
k,l−
1
2
,m
Φ
k,l,m+
1
2
−Φ
k,l,m−
1
2
'
(
)
)
)
)
*
+
,
,
,
,
=
1
H
∂Φ
k+
1
2
,l,m
∂x
−
∂Φ
k−
1
2
,l,m
∂x
+
∂Φ
k,l+
1
2
,m
∂y
−
∂Φ
k,l−
1
2
,m
∂y
+
∂Φ
k,l,m+
1
2
∂z
−
∂Φ
k,l,m−
1
2
∂z
'
(
)
)
*
+
,
,
=
1
H
2
Φ
k+1,l,m
−2Φ
k,l,m
+Φ
k−1,l,m
+Φ
k,l+1,m
−2Φ
k,l,m
+Φ
k,l−1,m
+Φ
k,l,m+1
−2Φ
k,l,m
+Φ
k,l,m−1
( )
obtain iterative solver by discretizing differential equation

46. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation relaxation technique
discretized Poisson’s equation
€
Φ
k,l,m
=
1
6
(Φ
k+1,l,m
+Φ
k−1,l,m
+Φ
k,l+1,m
+Φ
k,l−1,m
+Φ
k,l,m+1
+Φ
k,l,m−1
−ρ
k,l,m
H
2
)
€
ΔΦ
k,l,m
=ρ
k,l,m
Particle-Mesh Method
obtain iterative solver by discretizing differential equation

47. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation relaxation technique
discretized Poisson’s equation
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
€
ΔΦ
k,l,m
=ρ
k,l,m
Particle-Mesh Method
iterative solution:
€
Φ
k,l,m
i
→Φ
k,l,m
i+1
obtain iterative solver by discretizing differential equation

48. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
Particle-Mesh Method
!numerically integrate Poisson’s equation relaxation technique
k,l+1,m
k,l-1,m
k-1,l,mk+1,l,mk,l,m
obtain iterative solver by discretizing differential equation

49. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
Particle-Mesh Method
!numerically integrate Poisson’s equation relaxation technique
k,l+1,m
k,l-1,m
k-1,l,mk+1,l,mk,l,m
obtain iterative solver by discretizing differential equation
applicable to grids of
arbitrarygeometry!

50. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
Particle-Mesh Method
!numerically integrate Poisson’s equation relaxation technique
•how to sweep through the grid?
k,l+1,m
k,l-1,m
k-1,l,mk+1,l,mk,l,m

51. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
Particle-Mesh Method
!numerically integrate Poisson’s equation relaxation technique
•how to sweep through the grid?
?

52. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
Particle-Mesh Method
!numerically integrate Poisson’s equation relaxation technique
•how to sweep through the grid?
?

53. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
Particle-Mesh Method
k,l+1,m
k,l-1,m
k-1,l,mk+1,l,mk,l,m
!numerically integrate Poisson’s equation relaxation technique
•Gauss-Seidel sweeps:

54. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation relaxation technique
•Gauss-Seidel sweeps:
•loop over all “black” cells
•loop over all “red” cells
discretized Poisson’s equation
how many iterations iare necessary?
one iteration of the potential
€
Φ
k,l,m
i
→Φ
k,l,m
i+1
Particle-Mesh Method
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)

55. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:
discretized Poisson’s equation
Particle-Mesh Method
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
ΔΦk,l,mi=ρk,l,m

56. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:
discretized Poisson’s equation
Particle-Mesh Method
€
ΔΦ
k,l,m
i
−ρ
k,l,m
?
& → & 0
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)

57. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
€
ΔΦ
k,l,m
i
−ρ
k,l,m
?
& → & 0
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
density as given by mass assignment scheme!
density as given by currently best guess for !i !
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:

58. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
€
ΔΦ
k,l,m
i
−ρ
k,l,m
?
& → & 0
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)

€
•=suitable norm
€
R
i
=ΔΦ
k,l,m
i
−ρ
k,l,m
residual:
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:

59. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
€
ΔΦ
k,l,m
i
−ρ
k,l,m
?
& → & 0

€
•=suitable norm
error estimatetolerance
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
€
R
i
=ΔΦ
k,l,m
i
−ρ
k,l,m
≤εΤ
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:

60. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
€
R
i
=ΔΦ
k,l,m
i
−ρ
k,l,m
≤εΤ
!truncation error:
error due to discreteness of grid
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:

61. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
!truncation error:
error due to discreteness of grid
compare solution on actual grid
to solution on coarser grid

€
estimation
" → " " "
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
€
R
i
=ΔΦ
k,l,m
i
−ρ
k,l,m
≤εΤ
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:

62. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method

€
Τ
k,l,m
=PΔRΦ
k,l,m
i( )[ ]
−ΔΦ
k,l,m
i( )
€
RΦ
k,l,m
i
=Φ
j,n,p
i
ΔRΦk,l,mi( )=ρj,n,pi−1

€
PΔRΦ
k,l,m
i( )[ ]
=ρ
k,l,m
i
!truncation error:
restriction to coarser grid
prolongation to finer grid
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
Ri=ΔΦk,l,mi−ρk,l,m≤εΤ=εTk,l,m
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:

63. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
€
=ΔΦ
k,l,m
i
?
!truncation error:

€
PΔRΦ
k,l,m
i( )[ ]
=ρ
k,l,m
i
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:

€
Τ
k,l,m
=PΔRΦ
k,l,m
i( )[ ]
−ΔΦ
k,l,m
i( )
€
RΦ
k,l,m
i
=Φ
j,n,p
i
ΔRΦk,l,mi( )=ρj,n,pi−1
Ri=ΔΦk,l,mi−ρk,l,m≤εΤ=εTk,l,m

64. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
€
R
i
€
i
€
εT
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:
Ri=ΔΦk,l,mi−ρk,l,m≤εΤ=εTk,l,m

65. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
€
R
i
€
i
€
εT
faster convergence!?
!numerically integrate Poisson’s equation relaxation technique
•stopping criterion:
Ri=ΔΦk,l,mi−ρk,l,m≤εΤ=εTk,l,m

66. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation relaxation technique
•convergence:
discretized Poisson’s equation
Particle-Mesh Method
!slow convergence:
€
R
i+1
≈R
i
large-scale errors in "cannot be “relaxed”
sufficiently fast on the actual grid
=> use coarser grids to speed up convergence…
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
€
R
i
=ΔΦ
k,l,m
i
−ρ
k,l,m

67. Computational Astrophysics
Solving for Gravity
discretized Poisson’s equation
Particle-Mesh Method
!slow convergence:
€
R
i+1
≈R
i
multi-grid relaxation techniques
=> beyond the scope of this lecture though…
€
Φ
k,l,m
i+1
=
1
6
(Φ
k+1,l,m
i
+Φ
k−1,l,m
i
+Φ
k,l+1,m
i
+Φ
k,l−1,m
i
+Φ
k,l,m+1
i
+Φ
k,l,m−1
i
−ρ
k,l,m
H
2
)
€
R
i
=ΔΦ
k,l,m
i
−ρ
k,l,m
!numerically integrate Poisson’s equation relaxation technique
•convergence:

68. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
accuracy of either
relaxation
or
FFT method
to solve Poisson’s equation?
Particle-Mesh Method

69. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
pure PM calculation
Particle-Mesh Method

70. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
pure PM calculation
Particle-Mesh Method
distance measured in grid cells

71. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
pure PM calculation
Particle-Mesh Method
the force is automatically softened... (cf. tree-code lecture)
distance measured in grid cells

72. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
pure PM calculation
Particle-Mesh Method
the force is automatically softened...but what if need to resolve smaller scales?
distance measured in grid cells

73. Computational Astrophysics
Solving for Gravity
pure PM calculation
!numerically integrate Poisson’s equation
introduction of finer grids in high density regions
Particle-Mesh Method
distance measured in grid cells

74. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
Particle-Mesh Method
Yahagi & Yoshi (2001)
AMR calculation

75. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
Particle-Mesh Method
Yahagi & Yoshi (2001)
details later…
AMR calculation

76. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
pure PM calculation
Particle-Mesh Method
…and what are these wiggles?

77. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
pure PM calculation
Particle-Mesh Method
…and what are these wiggles?

78. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation
pure PM calculation
Particle-Mesh Method
force anisotropy
…and what are these wiggles?

79. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation

€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid

€
!
F (
!
g
k,l,m
)
3. differentiate potential to get forces

€
!
F (
!
g
k,l,m
)→
!
F (
!
x
i
)
4. interpolate forces back to particles
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)

80. Computational Astrophysics
Solving for Gravity
!obtaining the forces

€
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)

€
F
x
(
!
g
k,l,m
)=−m
Φ(
!
g
k+1,l,m
)−Φ(
!
g
k−1,l,m
)
2H
F
y
(
!
g
k,l,m
)=−m
Φ(
!
g
k,l+1,m
)−Φ(
!
g
k,l−1,m
)
2H
F
z
(
!
g
k,l,m
)=−m
Φ(
!
g
k,l,m+1
)−Φ(
!
g
k,l,m−1
)
2H

€
H=(current) grid spacing
Particle-Mesh Method
k,l+1,m
k,l-1,m
k-1,l,mk+1,l,mk,l,m

81. Computational Astrophysics
Solving for Gravity
!numerically integrate Poisson’s equation

€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid

€
!
F (
!
g
k,l,m
)
3. differentiate potential to get forces

€
!
F (
!
g
k,l,m
)→
!
F (
!
x
i
)
4. interpolate forces back to particles
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)
Particle-Mesh Method

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)

82. Computational Astrophysics
Solving for Gravity
!interpolating the forces

€
!
F (
!
g
k,l,m
)→
!
F (
!
r
i
)
Particle-Mesh Method

83. Computational Astrophysics
Solving for Gravity
!interpolating the forces
in practice the triple sum is “only” over 8 (CIC) or 27 (TSC) cells…
Particle-Mesh Method

€
!
F (
!
g
k,l,m
)→
!
F (
!
r
i
)
use the inverse of the mass assignment scheme
to insure momentum conservation and minimize force anisotropies
!F(!ri)= !F(!gk,l,m)W(!ri−!gk,l,m)m∑l∑k∑

84. Computational Astrophysics
Solving for Gravity
!interpolating the forces
Particle-Mesh Method
*check by calculating the total (periodic) force:
!Ftot=!F(!ri)i=1N∑ =i=1N∑ !F(!gk,l,m)W(!ri−!gk,l,m)m∑l∑k∑=...= mimjH3W(!ri−!gk,l,m)G(!gk,l,m−!g#k,#l,#m)W(!rj−!g#k,#l,#m)#k,#l,#m∑k,l,m∑j=1N∑i=1N∑=...=!0

€
!
F (
!
g
k,l,m
)→
!
F (
!
r
i
)
use the inverse of the mass assignment scheme
to insure momentum conservation and minimize force anisotropies
!F(!ri)= !F(!gk,l,m)W(!ri−!gk,l,m)m∑l∑k∑
Fx(!gk,l,m)=−mΦ(!gk+1,l,m)−Φ(!gk−1,l,m)2HΦ(!gk,l,m)= G#m∑#l∑#k∑ (!gk,l,m−!g#k,#l,#m)ρ(!g#k,#l,#m)ρ(!gk,l,m)=M(!gk,l,m)H3M(!gk,l,m)=miW(!ri−!gk,l,m)i=1N∑
anti-symmetric!
(because of invariance under coordinate inversion)
PM scheme:

85. Computational Astrophysics
Solving for Gravity
!Particle-Mesh (PM) method
Summary

€
Φ(
!
g
k,l,m
)
2. solve Poisson’s equation on grid

€
!
F (
!
g
k,l,m
)
3. differentiate potential to get forces

€
!
F (
!
g
k,l,m
)→
!
F (
!
x
i
)
4. interpolate forces back to particles
1. calculate mass density on grid

€
!
x
i
→ρ(
!
g
k,l,m
)

€
ΔΦ(
!
g
k,l,m
)=4πGρ(
!
g
k,l,m
)
!
F (
!
g
k,l,m
)=−m∇Φ(
!
g
k,l,m
)
anyone fancies to write a PM code as the project?

86. Computational Astrophysics
Solving for GravityAlexander Knebe, Universidad Autonoma de Madrid