Module on Computational Astrophysics

1. Module on Computational Astrophysics Jim Stone Department of Astrophysical Sciences 125 Peyton Hall : ph. 258-3815: jmstone@princeton.edu www.astro.princeton.edu/~jstone Lecture 1: Introduction to astrophysics, mathematics, and methods Lecture 2: Optimization, parallelization, modern methods Lecture 3: Particle-mesh methods Lecture 4: Particle-based hydro methods, future directions

2. Do we really need to compute force from every star for distant objects? Andromeda – 2 million light years away

3. Distance = 25 times size Solving the force problem with software -- tree codes

4. Organize particles into a tree. In Barnes-Hut algorithm, use a quadtree in 2D

5. In 3D, Barnes-Hut uses an octree

6. If angle subtended by the particles contained in any node of tree is smaller than some criterion, then treat all particles as one. Results in an Nlog(N) algorithm.

7. Alternative to Barnes-Hut is KD tree. • KD tree is binary - extremely efficient • Requires N to be power of 2 • Nnodes = 2N-1

8. Parallelizing tree code. Best strategy is to distribute particles across processors. That way, work of computing forces and integration is distributed across procs. Challenge is load balancing • Equal particles  equal work. • Solution: Assign costs to particles based on the work they do • Work unknown and changes with time-steps • Insight : System evolves slowly • Solution: Count work per particle, and use as cost for next time-step.

9. A Partitioning Approach: ORB • Orthogonal Recursive Bisection: • Recursively bisect space into subspaces with equal work • Work is associated with bodies, as before • Continue until one partition per processor • High overhead for large no. of processors

10. Another Approach: Costzones • Insight: Tree already contains an encoding of spatial locality. • Costzones is low-overhead and very easy to program

11. Methods so far… • Time integration • second-order leap frog • variable time steps • importance of keeping time symmetry • higher-order (Hermite) schemes • Force evaluation • direct summation (PP) • tree codes

12. Running a billion particles. Even with GRAPE, direct summation is limited to N ~ 106 Tree codes can reach 107-8 particles (with very efficient parallelization) Even more efficient methods are required to go beyond 109 particles. Particle-Mesh Methods

13. Particle-mesh methods. • Calculate the force on a particle from the gravitational potential F of all the particles, rather than using Newton’s Law of gravity to evaluate inter-particle forces directly • is a scalar function of x given by Poisson equation Where r is the mass density

14. Steps in a PM code. Assign particles to mesh to compute r Solve Poisson equation with appropriate BCs Compute force on particles Integrate motion Illustrate concepts with a 1D algorithm.

15. i-1 i i+1 1. Assigning particles to grid Discretize space into Ng zones in 1D Simplest method is to assign entire mass of particle to zone which contains it = Nearest Grid Point scheme Results in a very course distribution of ri.

16. r m/H H x Better method is particle-in-cell (PIC) Treat particles as having finite size, assume some distribution for the mass of the particle within this size. e.g., use top-hat distribution for particle mass: Then distribute mass of particles amongst zones according to overlap Can adopt more complex shapes for mass-distribution of particles, e.g. Gaussian, etc. i-1 i i+1

17. r0 ri-1riri+1rNg+1 F0 Fi-1FiFi+1FNg+1 2. Solving Poisson equation In 1D, Poisson equation is: elliptic PDE Discretize F on a mesh at cell-centers using Ng grid points Then centered FDE for Poisson equation is: For i=1, need F0; for i=Ng need FNg+1; -- Boundary Conditions

18. r0 ri-1riri+1rNg+1 F0 Fi-1FiFi+1FNg+1 Getting boundary conditions Periodic boundary conditions are most common: F0 = FNg; FNg+1 = F1 Otherwise, standard to use a multipole expansion. Often, the first term (monopole) is good enough (if boundaries are far from most of the mass)

19. The FDE for Poisson equation Can be written as a set of coupled linear equations, Ng eqns. Results in a tri-diagonal matrix. Very efficient algorithms are possible for “tri-di” systems (Hockney & Eastwood p185). Forward elimination/Back substitution allows system to be solved in only 4Ng operations.

20. r0 ri-1riri+1rNg+1 F0 Fi-1FiFi+1FNg+1 3. Computing Force on particles Once potential is know, force computed from Note forces are located at cell edges: Fi-3/2 Fi-1/2 Fi+1/2 Fi+3/2 Some sort of spatial interpolation is needed to get force at location of particles. Order of interpolation should be consistent with order of function used to distribute mass in Step 1.

21. Extension to 3D. For particle assignment, mass smoothing function must be sphere For boundary conditions, use spherical multipole expansion Solving Poisson equation: Discretize in 3D as

22. FDE in 3D becomes: Note this can again be written as a set of coupled linear equations NxNyNz equations in all. Very large sparse banded matrix of size NxNyNz byNxNyNz Even a very small problem (203) leads to large matrix 8000 x 8000 “Typical” size grid (1003) leads to 106 x 106 matrix Clearly need very efficient methods to solve matrix. Same issues apply to any method for solving elliptic PDEs.

23. Methods to solve matrix • Relaxation schemes. Guess solution, then “relax”. e.g. Multigrid • Sparse banded solvers, e.g. Conjugate Gradient methods • Fourier methods - very powerful, but require periodic BCs

24. Adaptive PM Using a fixed grid to solve the Poisson equation is not efficient in regions with structure on small length scales. Instead, adaptively embed fine grids in dense regions. Grid spacing in each fine level factor of 2 smaller than next coarsest level. Dynamic range is then 2Nlevels-1. For elliptic systems, solution on fine mesh alters solution on coarse mesh globally -- must solve on all levels simultaneously