Performance analysis of direct N-body algorithms for astrophysical simulations on distributed systems

Abstract

We discuss the performance of direct summation codes used in the simulation of astrophysical stellar systems on highly distributed architectures. These codes compute the gravitational interaction among stars in an exact way and have an O(N²) scaling with the number of particles. They can be applied to a variety of astrophysical problems, like the evolution of star clusters, the dynamics of black holes, the formation of planetary systems, and cosmological simulations. The simulation of realistic star clusters with sufficiently high accuracy cannot be performed on a single workstation but may be possible on parallel computers or grids. We have implemented two parallel schemes for a direct N-body code and we study their performance on general purpose parallel computers and large computational grids. We present the results of timing analyzes conducted on the different architectures and compare them with the predictions from theoretical models. We conclude that the simulation of star clusters with up to a million particles will be possible on large distributed computers in the next decade. Simulating entire galaxies however will in addition require new hybrid methods to speedup the calculation.

Introduction

Numerical methods for solving the classical astrophysical N-body problem have evolved in two main directions in recent years. On the one hand, approximated models like Fokker–Planck models [1], [2], gaseous models [3], and Monte Carlo models [4], [5], [6] have been applied to the simulation of globular clusters and galactic nuclei. These models permit to follow the global evolution of large systems along their lifetime but at the expense of moderate accuracy and resolution. The basic approach is to group particles according to their spatial distribution and use a truncated multipole expansion to evaluate the force exerted by the whole group instead of evaluating directly the contributions from the single particles. On the other hand, direct summation methods have been developed to accurately model the dynamics and evolution of collisional systems like dense star clusters. These codes compute all the inter-particle forces and are therefore the most accurate. They are also more general, as they can be used to simulate both low and high density regions. Their high accuracy is necessary when studying physical phenomena like mass segregation, core collapse, dynamical encounters, formation of binaries or higher order systems, ejection of high velocity stars, and runaway collisions. Direct methods have an O(N²) scaling with the number of stars and are therefore limited to smaller particle numbers compared to approximated methods. For this reason, so far they have only been applied to the simulation of moderate size star clusters. The simulation of globular clusters containing one million stars is still a challenge from the computational point of view, but it is an important astrophysical problem. It will provide insight in the complex dynamics of these collisional systems, in the microscopic and macroscopic physical phenomena, and it will help finding evidence of the controversial presence of a central black hole. The simulation of the evolution of these systems under the combined effects of gravity, stellar evolution, and hopefully hydrodynamics will allow to study the stellar population and to compare the results with observations.

The need to simulate ever larger systems and to include a realistic physical treatment of stars asks for a speedup in the most demanding part of the calculation: the gravitational dynamics. Significant improvement in the performance of direct codes can be obtained either by means of special purpose computers like GRAPE hardware [7] or of general purpose distributed systems [8]. In this work, we focus on the performance of direct N-body codes on distributed systems. Two main classes of algorithms can be used to parallelize direct summation N-body codes: replicated data and distributed data algorithms. In this work, we present a performance analysis of the two algorithms on different architectures: a Beowulf cluster, three supercomputers, and two computational grids. We provide theoretical predictions and actual measurements for the execution time on different platforms, allowing the choice of the best performing scheme for a given architecture.