Implimentation of MVS integrators to cosmological N

Ryuji Morishima (UCLA/JPL)
 N-body code: Gravity solver + Integrator
 Gravity solver must be fast and handle close encounters
 Special hardware (N2): Grape, GPS
 Tree (N log N): PKDGRAV, Gadget
 Integrator must take a large time step with good accuracy
 Bulirsch-Stoer
 Hermite: often used with Grape
 Mixed Variable Symplectic (MVS) integrators:
SYMBA, Mercury
 This work: Implementation of SYMBA to PKDGRAV
 Developed by Stadel (2001)
 Source is open in The astro-code wiki
 Tree gravity: 4th order multiple moments
 Adaptive to various parallel environments
(shared memory, mpi)
 Different functions and integrators
 Collisions (Richardson et al. 2000)
 SPH (Wadsley et al. 2004)
 Fragmentation (Leinhardt & Richardson 2005)
 SYMBA integrator (Morishima et al. 2010)
Up to 4th order (Hexadecapole)
Error estimation from
cosmological simulations
k-D Tree
Spatial binary tree
Spatial binary tree can reduce the higher order multi-pole moments
It is also as efficient as k-D tree in neighboring search.
Mater layer
• Controls overall flows of
Processor Set Tree (PST) layer
• Assigns tasks to processors
Parallel KD layers
• Executes tasks in each core
One needs to understand PST format
but not parallel primitives such as MPI
Machine dependent layers
• Interface to parallel
primitives (e.g. MPI)
Specialized for systems with a massive central body
Mixed variables: Cartesian and Keplerian co-ordinates
A large time step along Keplerian orbit
Time-reversible (no secular error)
 Handling close encounters:
 SYMBA (Duncan et al. 1998)
 Mercury (Chambers 1999)
 Most of N-body simulations for planet formation have been
performed by these two codes in last decade
 But both codes use N2 gravity calculations
 Democratic co-ordinate
 Heliocentric position + barycentric velocity
Hkep>>Hint (if there is no close encounter)
 Potential (or Force) decomposition based on mutual
distance normalized by the Hill radius
 A higher order potential component is calculated with
a small block-sized time step
Kepler Drift with F1
Kick (F0)
Kick (F0)
 The time step size needs to be determined by the
minimum mutual distance during particle drift.
 This distance must be estimated by using particle coordinates at the beginning and ending of particle drift
symmetrically (e.g. Hut et al. 1995).
Half kick (t0/2) due to Sun’s motion
Half Kick (t0/2) due to force F0 from other particles
Kepler drift (t0) for all particles
Tree build and neighboring search (after drift)
Particles in close encounters
Sent back to pre-drift positions and velocities
Put into a single core (domain decomposition)
SYMBA multiple time stepping
Collisions are also handled here
Tree build and gravity calculation and neighboring
search (before drift)
7. 2 and 1
Martian meteorites
Accretion truncated
at 14 My
 In 3-body encounter, time stepping of these 3 bodies is
synchronized (for time symmetry)
 If the system’s number density is high, all particles
share the time step….
 PKDGRAV-SYMBA works as desired unless the
system’s number density is so high that most of bodies
are in multi-body encounters.
 For such systems, time symmetry needs to be

similar documents