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 program call Processor Set Tree (PST) layer • Assigns tasks to processors call 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) Hkep>>Hsun 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). 1. 2. 3. 4. 5. 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 1. 2. 3. 4. 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 6. With gas Without gas Martian meteorites Accretion truncated at 14 My Accretion Extrapolated 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 sacrificed?