Module on Computational Astrophysics

1
Module on Computational Astrophysics

• Jim Stone
• Department of Astrophysical Sciences
• 125 Peyton Hall ph. 258-3815
  jmstone_at_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
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 Newtons 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.

1. Assign particles to mesh to compute r
2. Solve Poisson equation with appropriate BCs
3. Compute force on particles
4. Integrate motion

Illustrate concepts with a 1D algorithm.
15
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
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 i1
17
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
i1, need F0 for iNg need FNg1 -- Boundary
Conditions
18
Getting boundary conditions
Periodic boundary conditions are most common
F0 FNg FNg1 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
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 Fi1/2 Fi3/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.

1. For particle assignment, mass smoothing function
   must be sphere
2. For boundary conditions, use spherical multipole
   expansion
3. 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
by NxNyNz 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
25
(No Transcript)
26
(No Transcript)
27
(No Transcript)