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

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

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

Motivating example from computational astrophysics 3

Example: Stellar dynamics in a globular cluster 4

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

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

Systems innovation! Who needs computers? 7

8

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

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

Example: Dark matter filament formation Particles = dark matter Dim: 43 Mpc ≈ 10 21 km Source: Andrey Kravtsov, U. Chicago – http://cosmicweb.uchicago.edu 11

Length-/time-scale challenges Long length-scales Star cluster is 10 9 times larger than star Neutron star: 10 14 x Long time-scales Close passage between stars ~ hours Evolution of cluster ~ 10 9 years ⇒ 10 14 x (neutron stars: 10 20 ) 12

Newton’s laws of motion and gravity d 2 r i ( t ) = m i F i dt 2 m i ¨ ≡ r i r i − r j � = − Gm i m j | r i − r j | 3 j � = i ⇓ r i − r j � ¨ = r i − G m j | r i − r j | 3 j � = i 13

Rewrite as 1st-order ODE F i r i − r j � − G m j = ¨ r i = | r i − r j | 3 m i j � = i ⇓ � � � � d r i v i = F i /m i v i dt n particles ⇒ 6 n -element vector 14

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(n 2 ) → O(n log n) → O(n) algorithms Applies to other domains, e.g. , molecular dynamics, … 15

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

ODE solvers 17

Recall: Computational tasks Setup Solve initial-value problem (ODE) Evaluate forces Integrate over some time step 18

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

Solve initial-value problem (ODE), with force evaluation at each step � � � � d r i v i = F i /m i v i dt ⇓ F i r i − r j � − G m j = ¨ r i = m i | r i − r j | 3 j � = i 20

“Soften” to remove singularities � � � � d r i v i = F i /m i v i dt ⇓ F i r i − r j � = ¨ = − G m j r i 3 m i ( | r i − r j | 2 + ǫ 2 ) 2 j � = i 21

IVP ODE solvers Many algorithms Design space Taylor series ( e.g. , Euler) Single vs. multistep Runge-Kutta Explicit vs. implicit Extrapolation No. of req’d function/ derivative evaluations Multistep Multivalue 22

IVP ODE solvers Many algorithms Design space Single vs. multistep Taylor series ( e.g. , Euler) Explicit vs. implicit Runge-Kutta No. of req’d function/ Extrapolation derivative evaluations Multistep Multivalue 23

Taylor series-based methods Derive finite-difference approximations using Taylor series ODE: ˙ = f ( t, y ( t )) y y ( t ) + h 2 y ( t + h ) = y ( t ) + h ˙ 2 ¨ y ( t ) + · · · Example: Forward Euler = t 0 + h · k k ∈ { 1 , 2 , . . . } t k y k + h · f ( t k , y k ) y k +1 ← 24

Taylor series-based methods: Forward Euler = t 0 + h · k k ∈ { 1 , 2 , . . . } t k y k + h · f ( t k , y k ) y k +1 ← 1st-order accurate (Higher order: Take more derivative terms) Single-step (variably-sized steps possible) Explicit : y k+1 depends only on t k Narrow stability window, i.e. , small range of values for h for stability 25

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

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

Taylor series-based methods: Backward Euler y ( t ) + h 2 y ( t − h ) = y ( t ) − h ˙ 2 ¨ y ( t ) + · · · y k − 1 + h · f ( t k , y k ) y k ← Implicit : y k+1 depends only on t k Larger stability window vs. forward Euler 1st-order accurate (Higher order: Take more derivative terms) Single-step (variably-sized steps possible) 28

IVP ODE solvers Many algorithms Design space Taylor series ( e.g. , Euler) Single vs. multistep Explicit vs. implicit Runge-Kutta No. of req’d function/ Extrapolation derivative evaluations Multistep Multivalue 29

Runge-Kutta methods Single-step, but w/ high-order derivs replaced by FD approx in [ t k , t k+1 ] Classical 4th-order RK: y k + h k 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) y k +1 ← = f ( t k , y k ) k 1 f ( t k + h k 2 , y k + h k = 2 k 1 ) k 2 f ( t k + h k 2 , y k + h k = 2 k 2 ) k 3 = f ( t k + h k , y k + h k k 3 ) k 4 30

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

IVP ODE solvers Many algorithms Design space Taylor series ( e.g. , Euler) Single vs. multistep Runge-Kutta Explicit vs. implicit No. of req’d function/ Extrapolation derivative evaluations Multistep Multivalue 32

IVP ODE solvers Many algorithms Design space Taylor series ( e.g. , Euler) Single vs. multistep Runge-Kutta Explicit vs. implicit Extrapolation No. of req’d function/ derivative evaluations Multistep Multivalue 33

Multistep methods “Memory:” Use information at more than one prior time step General form of linear multistep methods: m m � � α i y k +1 − i + h β i f ( t k +1 − i , y k +1 − i ) y k +1 ← i =1 i =0 Obtain α , β coefficients via polynomial interpolation β 0 = 0 ⇒ explicit; otherwise, implicit Predictor-corrector: Explicit guesses, implicit improvements 34

Multistep method examples Simplest 2nd-order explicit two-step method y k + h 2 (3 ˙ y k − ˙ y k − 1 ) y k +1 ← Simplest 2nd-order implicit two-step method y k + h = 2 ( ˙ y k +1 + ˙ y k ) y k +1 4th-order Adams-Bashforth predictor y k + h 24(55 ˙ y k − 59 ˙ y k − 1 − 37 ˙ y k − 2 − 9 ˙ y k − 3 ) y k +1 ← 4th-order implicit Adams-Moulton corrector y k + h = 24(9 ˙ y k +1 + 19 ˙ y k − 5 ˙ y k − 1 + ˙ y k − 2 ) y k +1 35

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

IVP ODE solvers Many algorithms Design space Taylor series ( e.g. , Euler) Single vs. multistep Runge-Kutta Explicit vs. implicit Extrapolation No. of req’d function/ derivative evaluations Multistep Multivalue 37

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

IVP ODE solvers Many algorithms Design space Taylor series ( e.g. , Euler) Single vs. multistep Runge-Kutta Explicit vs. implicit Extrapolation No. of req’d function/ derivative evaluations Multistep Multivalue 39

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

Administrivia 41

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

Parallel ODEs: Waveform relaxation 43

Picard iteration Consider IVP , a system of n-ODEs = f ( t, y ) y ˙ ≥ t t 0 y ( t 0 ) = y 0 Idea: Following iteration converges to soln. of IVP (if f is Lipschitz) � t Y k +1 ( t ) y 0 + f ( s, Y k ( s )) ds ← t 0 n independent 1-D quadrature calculations! 44

Jacobi waveform relaxation Decouple system into n independent ODEs Consider n = 2   � � t, y ( k +1) , y ( k )    y ( k +1) f 1 1 2 d 1   ←  dt y ( k +1) � �   t, y ( k ) 1 , y ( k +1) f 2 2 2 Red-black Gauss-Seidel (and other) splittings possible 45

Back to the n-body problem 46

Recall: Computational tasks Setup Solve IVP Evaluate forces Integrate over some time step 47

Leap-frog integrator ¨ r ( t ) = g ( r ( t )) Consider discretized position and velocity on a staggered grid: r 0 r 1 r 2 r 3 r 4 ˙ ˙ ˙ ˙ ˙ r 1 r 3 r 5 r 7 r 9 2 2 2 2 2 h 48