1 / 33

Module on Computational Astrophysics

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

lucius-farmer
Download Presentation

Module on Computational Astrophysics

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  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. Computational Tasks Setup distribution of N particles Compute forces between particles Evolve positions using ODE solver Display/analyze results

  3. n-1 n-1/2 n n+1/2 n+1 time x, F v x, F v x, F Leap-frog scheme Define positions (x) and forces (F) at time level n velocities (v) at time level n+1/2 Then, for ith particle To start integration, need initial x and V at two separate time levels. Specify x0 and v0 and then integrate V to Dt/2 using high-order scheme

  4. Accuracy of leap-frog scheme Can show the truncation error in leap-frog is second order in Dt Evolution eqns: Replace Vn-1/2 in second equation using first Substitute this back into first equation Rearrange This is central difference formula for F=ma.

  5. Let X(t) be the “true” (analytic) solution. Then EQ. 1 Use a Taylor expansion to compute Xn+1 and Xn-1 Thus Substitute these back into eq. 1 Truncation error O(Dt2)

  6. Truncation versus round-off error Note the error we have just derived is truncation error Unavoidable result of approximating solution to some order in x or t Completely unrelated to round-off error, which results from representing the continuous set of real numbers with a finite number of bits. Truncation error can be reduced by using smaller step Dt, or higher-order algorithm Round-off error can be reduced by using higher precision (64 bit rather than 32, etc.), and by ordering operations carefully. In general, truncation error is much larger than round-off error

  7. Stability of leap-frog scheme Easiest to illustrate with an example. Suppose the force is given by a harmonic oscillator, that is: Then “true” (analytic) solution is Substitute force law into leap-frog FDE for F=ma Look for oscillatory solutions of the form x = x0eiwt giving

  8. This is good! Leap-frog (correctly) gets oscillatory solutions, but at a modified frequency Note this gives correct solution (W) as Dt --> 0 For WDt > 2, frequency becomes complex. Real part of w’ gives oscillatory solution, imaginary part gives exponentially growing (unstable) solution. So stability limit is Dt < 2/W; or Dt < 2/[(dF/dx)/m]1/2 in general Above is a simple example of a von-Neumann stability analysis

  9. n-1 n-1/2 n n+1/2 n+1 time x, F v x, F v x, F Consistency of leap-frog scheme Leap-frog is consistent in the sense that as Dt --> 0, the difference equations converge to the differential equations Leap-frog is also a symplectic method (time symmetric). Scheme has same accuracy for Dt negative.

  10. Efficiency of leap-frog scheme Leap-frog is extremely efficient in terms of computational cost (only 12 flops per particle excluding force evaluation) Also extremely efficient in terms of memory storage (does not require storing multiple time levels). All the work (and memory) is in force evaluation: 10N flops per particle for direct summation To update all particle positions in one second on a 1 Gflop processor requires N < 104 Extra efficiency can be gained by using different timesteps for each particle (more later).

  11. 4. Display and Analysis Display requires plotting particles as spatial points in 3D. Sophisticated packages render different particles with different colors, allow perspective to be rotated, fly-throughs, etc.

  12. Quantitative analysis requires following diagnostics such as E, L, etc. for all particles. Need to store 6 words of data for each particle (x and v), in enough data dumps to give good time resolution of dynamics. Biggest simulations follow 109 particles, 1000 data dumps requires 48 TB (64 bit words).

  13. Putting it all together. Should be able to put together a direct N-body code using leap-frog. • Initialize x, V for N particles in region of size R • Use constant time step which is << crossing time • Use direct summation to compute F, with softened potential • Update x,V with leap-frog for a few crossing times • Output x for each particle frequently compared to crossing time • Use gnuplot, etc. to plot motion of particles • Experiment with different N

  14. Or, you can just download the code from the web… Large collection of N-body codes available from Peter Teuben at University of Maryland http://bima.astro.umd.edu/nemo This include Sverre Aarseth’s original NBODY0 code (first widely used direct N-body (PP) code used in astrophysics) Also see http://www.artcompsci.org/vol_1/v1_web/v1_web.html

  15. /* * LEAPINT.C: program to integrate hamiltonian system using leapfrog. */ #include <math.h> #define MAXPNT 100 /* maximum number of points */ void leapstep(); /* routine to take one step */ void accel(); /* accel. for harmonic osc. */ void printstate(); /* print out system state */ void main(argc, argv) int argc; char *argv[]; { int n, mstep, nout, nstep; double x[MAXPNT], v[MAXPNT], tnow, dt; /* first, set up initial conditions */ n = 1; /* set number of points */ x[0] = 1.0; /* set initial position */ v[0] = 0.0; /* set initial velocity */ tnow = 0.0; /* set initial time */ /* next, set integration parameters */ mstep = 256; /* number of steps to take */ nout = 4; /* steps between outputs */ dt = 1.0 / 32.0; /* timestep for integration */ /* now, loop performing integration */ for (nstep = 0; nstep < mstep; nstep++) { /* loop mstep times in all */ if (nstep % nout == 0) /* if time to output state */ printstate(x, v, n, tnow); /* then call output routine */ leapstep(x, v, n, dt); /* take integration step */ tnow = tnow + dt; /* and update value of time */ } if (mstep % nout == 0) /* if last output wanted */ printstate(x, v, n, tnow); /* then output last step */ }

  16. /* * LEAPSTEP: take one step using the leap-from integrator, formulated * as a mapping from t to t + dt. WARNING: this integrator is not * accurate unless the timestep dt is fixed from one call to another. */ void leapstep(x, v, n, dt) double x[]; /* positions of all points */ double v[]; /* velocities of all points */ int n; /* number of points */ double dt; /* timestep for integration */ { int i; double a[MAXPNT]; accel(a, x, n); /* call acceleration code */ for (i = 0; i < n; i++) /* loop over all points... */ v[i] = v[i] + 0.5 * dt * a[i]; /* advance vel by half-step */ for (i = 0; i < n; i++) /* loop over points again...*/ x[i] = x[i] + dt * v[i]; /* advance pos by full-step */ accel(a, x, n); /* call acceleration code */ for (i = 0; i < n; i++) /* loop over all points... */ v[i] = v[i] + 0.5 * dt * a[i]; /* and complete vel. step */ }

  17. /* * ACCEL: compute accelerations for harmonic oscillator(s). */ void accel(a, x, n) double a[]; /* accelerations of points */ double x[]; /* positions of points */ int n; /* number of points */ { int i; for (i = 0; i < n; i++) /* loop over all points... */ a[i] = - x[i]; /* use linear force law */ } /* * PRINTSTATE: output system state variables. */ void printstate(x, v, n, tnow) double x[]; /* positions of all points */ double v[]; /* velocities of all points */ int n; /* number of points */ double tnow; /* current value of time */ { int i; for (i = 0; i < n; i++) /* loop over all points... */ printf("%8.4f%4d%12.6f%12.6f\n", tnow, i, x[i], v[i]); }

  18. Variable time steps with leap-frog For efficiency, need to take variable time steps (evolve particles at center of cluster on smaller timestep than particles at edge).

  19. n-1 n-1/2 n n+1/2 n+1 time x, F v x, F v x, F Variable time steps with leap-frog However, this destroys symmetry of leap-frog; greatly increases truncation error.

  20. But, variable timestep leap-frog can be symmetrized Hut, Makino, & McMillan 1995

  21. Force evaluation with variable time steps. Now particle positions are known at different time levels. Greatly complicates force calculation. Must compute derivatives of force wrt time, and use Taylor expansion to compute total force on particle at current position. The Good: Allows higher-order (Hermite) integration methods. The Bad: This just makes force evaluation even more expensive! The Ugly: Direct N-body must be optimized if we are to go beyond 104 particles.

  22. Jun Makino, U. Tokyo Solving the force problem with hardware. Special purpose hardware to compute force:

  23. GRAPE-6 • The 6th generation of GRAPE (Gravity Pipe) Project • Gravity calculation with 31 Gflops/chip • 32 chips / board ⇒ 0.99 Tflops/board • 64 boards of full system is installed in University of Tokyo ⇒ 63 Tflops • On each board, all particle data are set onto SRAM memory, and each target particle data is injected into the pipeline, then acceleration data is calculated • Gordon Bell Prize at SC2000, SC2001(Prof. Makino, U. Tokyo)also nominated at SC2002

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

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

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

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

  28. 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.

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

  30. 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.

  31. 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

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

  33. Space Filling Curves Peano-Hilbert Order Morton Order

More Related