[PPT] - The n-body problem (1/2) Prof. Richard Vuduc Georgia Institute of PowerPoint Presentation

SLIDE 1

The n-body problem (1/2)

Prof. Richard Vuduc

Georgia Institute of Technology CSE/CS 8803 PNA: Parallel Numerical Algorithms [L.20] Tuesday, March 25, 2008

1

SLIDE 2

Today’s sources

CS 267 at UCB (Demmel & Yelick) Computational Science and Engineering, by Gilbert Strang Lectures 18–21 of M. Ricotti’s Astro 415 course at U. Maryland Jim Stone (Princeton) Andrey Kravtsov (U. Chicago) Mike Heath (UIUC) Changa (ChaNGa, U. Washington; based on CHARM++)

2

SLIDE 3

Motivating example from computational astrophysics

3

SLIDE 4

Example: Stellar dynamics in a globular cluster

4

SLIDE 5

Source: Frank Summers, AMNH, via Andrey Kravtsov Example: Stellar dynamics in a globular cluster

5

SLIDE 6

A little history of the n-body problem...

6

SLIDE 7

Systems innovation! Who needs computers?

7

SLIDE 8

8

SLIDE 9

Source: Ben Moore, U. Zurich (via Andrey Kravtsov) – http://nbody.net Example: Formation of group of galaxies 100 million light years ≈ 1020 km

9

SLIDE 10

Source: Andrey Kravtsov, U. Chicago – http://cosmicweb.uchicago.edu Example: Formation of group of galaxies

10

SLIDE 11

Source: Andrey Kravtsov, U. Chicago – http://cosmicweb.uchicago.edu Example: Dark matter filament formation

Particles = dark matter Dim: 43 Mpc ≈1021 km

11

SLIDE 12

Length-/time-scale challenges

Long length-scales

Star cluster is 109 times larger than star Neutron star: 1014 x

Long time-scales

Close passage between stars ~ hours Evolution of cluster ~ 109 years ⇒ 1014 x (neutron stars: 1020)

12

SLIDE 13

Newton’s laws of motion and gravity

Fi = mi d2ri(t) dt2 ≡ mi¨ ri = −Gmi

j=i

mj ri − rj |ri − rj|3 ⇓ ¨ ri = −G

j=i

mj ri − rj |ri − rj|3

13

SLIDE 14

Rewrite as 1st-order ODE

n particles ⇒ 6n-element vector

Fi mi = ¨ ri = −G

j=i

mj ri − rj |ri − rj|3 ⇓ d dt

ri

vi

=
vi

Fi/mi

14

SLIDE 15

Computational astrophysics

Sample problems

Planetary orbits over billions of years Planetary, star, and galaxy formation Star cluster evolution over time

Core computation: N-body problem

Integrate ODE – Off-the-shelf methods Inter-particle force evaluation – O(n2) → O(n log n) → O(n) algorithms

Applies to other domains, e.g., molecular dynamics, …

15

SLIDE 16

Source: James Stone (Princeton) Progress in N-body simulations for globular clusters Doubling time ≈ 3 years

16

SLIDE 17

ODE solvers

17

SLIDE 18

Recall: Computational tasks

Setup Solve initial-value problem (ODE)

Evaluate forces Integrate over some time step

18

SLIDE 19

Source: Jim Stone (Princeton) Setup initial distribution of particles

19

SLIDE 20

d dt

ri

vi

=
vi

Fi/mi

⇓

Fi mi = ¨ ri = −G

j=i

mj ri − rj |ri − rj|3

Solve initial-value problem (ODE), with force evaluation at each step

20

SLIDE 21

“Soften” to remove singularities

d dt

ri

vi

=
vi

Fi/mi

⇓

Fi mi = ¨ ri = −G

j=i

mj ri − rj (|ri − rj|2 + ǫ2)

3 2

21

SLIDE 22

IVP ODE solvers

Many algorithms

Taylor series (e.g., Euler) Runge-Kutta Extrapolation Multistep Multivalue

Design space

Single vs. multistep Explicit vs. implicit

No. of req’d function/

derivative evaluations

22

SLIDE 23

IVP ODE solvers

Many algorithms

Taylor series (e.g., Euler) Runge-Kutta Extrapolation Multistep Multivalue

Design space

Single vs. multistep Explicit vs. implicit

No. of req’d function/

derivative evaluations

23

SLIDE 24

Taylor series-based methods

Derive finite-difference approximations using Taylor series Example: Forward Euler

ODE: ˙ y = f (t, y(t)) y(t + h) = y(t) + h ˙ y(t) + h2 2 ¨ y(t) + · · · tk = t0 + h · k k ∈ {1, 2, . . .} yk+1 ← yk + h · f(tk, yk)

24

SLIDE 25

Taylor series-based methods: Forward Euler

1st-order accurate (Higher order: Take more derivative terms) Single-step (variably-sized steps possible) Explicit: yk+1 depends only on tk Narrow stability window, i.e., small range of values for h for stability

tk = t0 + h · k k ∈ {1, 2, . . .} yk+1 ← yk + h · f(tk, yk)

25

SLIDE 26

Forward Euler: Local and global errors Source: M. Heath (UIUC)

26

SLIDE 27

Forward Euler: Local and global errors Source: M. Heath (UIUC)

27

SLIDE 28

Taylor series-based methods: Backward Euler

Implicit: yk+1 depends only on tk Larger stability window vs. forward Euler 1st-order accurate (Higher order: Take more derivative terms) Single-step (variably-sized steps possible)

y(t − h) = y(t) − h ˙ y(t) + h2 2 ¨ y(t) + · · · yk ← yk−1 + h · f(tk, yk)

28

SLIDE 29

IVP ODE solvers

Many algorithms

Taylor series (e.g., Euler) Runge-Kutta Extrapolation Multistep Multivalue

Design space

Single vs. multistep Explicit vs. implicit

No. of req’d function/

derivative evaluations

29

SLIDE 30

Runge-Kutta methods

Single-step, but w/ high-order derivs replaced by FD approx in [tk, tk+1] Classical 4th-order RK:

yk+1 ← yk + hk 6 (k1 + 2k2 + 2k3 + k4) k1 = f(tk, yk) k2 = f(tk + hk 2 , yk + hk 2 k1) k3 = f(tk + hk 2 , yk + hk 2 k2) k4 = f(tk + hk, yk + hkk3)

30

SLIDE 31

Runge-Kutta methods

Nice properties

Self-starting Easy to program

Possible downsides

No a prior guidance on step-size selection Several function evaluations

“Workarounds:” Embedded RK methods; automation; implicit versions

31

SLIDE 32

IVP ODE solvers

Many algorithms

Taylor series (e.g., Euler) Runge-Kutta Extrapolation Multistep Multivalue

Design space

Single vs. multistep Explicit vs. implicit

No. of req’d function/

derivative evaluations

32

SLIDE 33

IVP ODE solvers

Many algorithms

Taylor series (e.g., Euler) Runge-Kutta Extrapolation Multistep Multivalue

Design space

Single vs. multistep Explicit vs. implicit

No. of req’d function/

derivative evaluations

33

SLIDE 34

Multistep methods

“Memory:” Use information at more than one prior time step General form of linear multistep methods: Obtain α, β coefficients via polynomial interpolation β0 = 0 ⇒ explicit; otherwise, implicit Predictor-corrector: Explicit guesses, implicit improvements

yk+1 ←

m

i=1

αiyk+1−i + h

m

i=0

βif(tk+1−i, yk+1−i)

34

SLIDE 35

Multistep method examples

Simplest 2nd-order explicit two-step method Simplest 2nd-order implicit two-step method 4th-order Adams-Bashforth predictor 4th-order implicit Adams-Moulton corrector

yk+1 ← yk + h 2 (3 ˙ yk − ˙ yk−1) yk+1 ← yk + h 24(55 ˙ yk − 59 ˙ yk−1 − 37 ˙ yk−2 − 9 ˙ yk−3) yk+1 = yk + h 24(9 ˙ yk+1 + 19 ˙ yk − 5 ˙ yk−1 + ˙ yk−2) yk+1 = yk + h 2 ( ˙ yk+1 + ˙ yk)

35

SLIDE 36

Multistep methods tradeoffs

Upsides

Good local error estimates via (predictor) - (corrector) Interpolation-based ⇒ results at points other than integration points Well-designed implicit multistep good for “stiff” problems

Downsides

Not self-starting Complicated to change step size Complicated to program

36

SLIDE 37

IVP ODE solvers

Many algorithms

Taylor series (e.g., Euler) Runge-Kutta Extrapolation Multistep Multivalue

Design space

Single vs. multistep Explicit vs. implicit

No. of req’d function/

derivative evaluations

37

SLIDE 38

Multivalue methods

Multistep with variable (adaptive) time steps Example: 4-value method y(t + h) = y(t) + h ˙ y(t) + h2 2 ¨ y(t) + h3 6 ... y (t) + . . . Yk ≡     yk h ˙ yk h2/2 · ¨ yk h3/6 · ... y k     Yk+1 ←     1 1 1 1 1 2 3 1 3 1     · Yk + α

38

SLIDE 39

IVP ODE solvers

Many algorithms

Taylor series (e.g., Euler) Runge-Kutta Extrapolation Multistep Multivalue

Design space

Single vs. multistep Explicit vs. implicit

No. of req’d function/

derivative evaluations

39

SLIDE 40

Sources of parallelism

Multistage, e.g., Runge-Kutta: Stage evaluation Evaluating right-hand side, e.g., forces in gravitational n-body Solving linear/non-linear systems, e.g., implicit methods Partitioning equations in multiple tasks – next

40

SLIDE 41

Administrivia

41

SLIDE 42

Upcoming schedule changes

Some adjustment of topics (TBD) Tu 4/1 — No class Th 4/3 — Attend talk by Doug Post from DoD HPC Modernization Program

42

SLIDE 43

Parallel ODEs: Waveform relaxation

43

SLIDE 44

Picard iteration

Consider IVP , a system of n-ODEs Idea: Following iteration converges to soln. of IVP (if f is Lipschitz) n independent 1-D quadrature calculations!

˙ y = f(t, y) t ≥ t0 y(t0) = y0 Yk+1(t) ← y0 + t

t0

f(s, Yk(s))ds

44

SLIDE 45

Jacobi waveform relaxation

Decouple system into n independent ODEs Consider n = 2 Red-black Gauss-Seidel (and other) splittings possible

d dt   y(k+1)

1

y(k+1)

2

  ←    f1

t, y(k+1)

1

, y(k)

2

f2
t, y(k)

1 , y(k+1) 2



 

45

SLIDE 46

Back to the n-body problem

46

SLIDE 47

Recall: Computational tasks

Setup Solve IVP

Evaluate forces Integrate over some time step

47

SLIDE 48

Leap-frog integrator

Consider discretized position and velocity on a staggered grid:

h ˙ r 3

2

˙ r 1

2

˙ r 5

2

˙ r 7

2

˙ r 9

2

r4 r3 r2 r1 r0

¨ r(t) = g(r(t))

48

SLIDE 49

Leap-frog integrator

vk+ 1

2 ≡ ˙

r

t + h

2

rk+1

← rk + h · vk+ 1

2

vk+ 3

2

← vk+ 1

2 + h · g(rk+1)

h ˙ r 3

2

˙ r 1

2

˙ r 5

2

˙ r 7

2

˙ r 9

2

r4 r3 r2 r1 r0

49

SLIDE 50

Time reversibility of leap-frog

rk+1 ← rk + h · vk+ 1

2

vk+ 3

2

← vk+ 1

2 + h · g(rk+1)

⇓ vk+ 1

2

← vk+ 3

2 − h · g(rk+1)

rk ← rk+1 − h · vk+ 1

2

50

SLIDE 51

Time reversibility of leap-frog Source: http://www.physics.drexel.edu/courses/Comp_Phys/Integrators/leapfrog/

51

SLIDE 52

“In conclusion…”

52

SLIDE 53

Backup slides

53