GADGET: a code for collisionless and gasdynamical cosmological simulations

Abstract

We describe the newly written code GADGET which is suitable both for cosmological simulations of structure formation and for the simulation of interacting galaxies. GADGET evolves self-gravitating collisionless fluids with the traditional N-body approach, and a collisional gas by smoothed particle hydrodynamics. Along with the serial version of the code, we discuss a parallel version that has been designed to run on massively parallel supercomputers with distributed memory. While both versions use a tree algorithm to compute gravitational forces, the serial version of GADGET can optionally employ the special-purpose hardware GRAPE instead of the tree. Periodic boundary conditions are supported by means of an Ewald summation technique. The code uses individual and adaptive timesteps for all particles, and it combines this with a scheme for dynamic tree updates. Due to its Lagrangian nature, GADGET thus allows a very large dynamic range to be bridged, both in space and time. So far, GADGET has been successfully used to run simulations with up to 7.5×10⁷ particles, including cosmological studies of large-scale structure formation, high-resolution simulations of the formation of clusters of galaxies, as well as workstation-sized problems of interacting galaxies. In this study, we detail the numerical algorithms employed, and show various tests of the code. We publicly release both the serial and the massively parallel version of the code.

Introduction

Numerical simulations of three-dimensional self-gravitating fluids have become an indispensable tool in cosmology. They are now routinely used to study the non-linear gravitational clustering of dark matter, the formation of clusters of galaxies, the interactions of isolated galaxies, and the evolution of the intergalactic gas. Without numerical techniques the immense progress made in these fields would have been nearly impossible, since analytic calculations are often restricted to idealized problems of high symmetry, or to approximate treatments of inherently nonlinear problems.

The advances in numerical simulations have become possible both by the rapid growth of computer performance and by the implementation of ever more sophisticated numerical algorithms. The development of powerful simulation codes still remains a primary task if one wants to take full advantage of new computer technologies.

Early simulations (Aarseth et al., 1979, Holmberg, 1941, Peebles, 1970, Press and Schechter, 1974, White, 1976, White, 1978; among others) largely employed the direct summation method for the gravitational N-body problem, which remains useful in collisional stellar dynamical systems, but it is inefficient for large N due to the rapid increase of its computational cost with N. A large number of groups have therefore developed N-body codes for collisionless dynamics that compute the large-scale gravitational field by means of Fourier techniques. These are the PM, P³M, and AP³M codes (Bertschinger and Gelb, 1991, Couchman, 1991, Eastwood and Hockney, 1974, Efstathiou et al., 1985, Hockney and Eastwood, 1981, Hohl, 1978, MacFarland et al., 1998). Modern versions of these codes supplement the force computation on scales below the mesh size with a direct summation, and/or they place mesh refinements on highly clustered regions. Poisson’s equation can also be solved on a hierarchically refined mesh by means of finite-difference relaxation methods, an approach taken in the ART code by Kravtsov et al. (1997).

An alternative to these schemes are the so-called tree algorithms, pioneered by Appel, 1981, Appel, 1985. Tree algorithms arrange particles in a hierarchy of groups, and compute the gravitational field at a given point by summing over multipole expansions of these groups. In this way the computational cost of a complete force evaluation can be reduced to a $\text{O}\text{(N}\text{log}\text{N)}$ scaling. The grouping itself can be achieved in various ways, for example with Eulerian subdivisions of space (Barnes and Hut, 1986), or with nearest-neighbour pairings (Jernigan and Porter, 1989, Press, 1986). A technique related to ordinary tree algorithms is the fast multipole-method (e.g., Greengard and Rokhlin, 1987), where multipole expansions are carried out for the gravitational field in a region of space.

While mesh-based codes are generally much faster for close-to-homogeneous particle distributions, tree codes can adapt flexibly to any clustering state without significant losses in speed. This Lagrangian nature is a great advantage if a large dynamic range in density needs to be covered. Here tree codes can outperform mesh based algorithms. In addition, tree codes are basically free from any geometrical restrictions, and they can be easily combined with integration schemes that advance particles on individual timesteps.

Recently, PM and tree solvers have been combined into hybrid Tree-PM codes (Bagla, 1999, Bode et al., 2000, Xu, 1995). In this approach, the speed and accuracy of the PM method for the long-range part of the gravitational force are combined with a tree-computation of the short-range force. This may be seen as a replacement of the direct summation PP part in P³M codes with a tree algorithm. The Tree-PM technique is clearly a promising new method, especially if large cosmological volumes with strong clustering on small scales are studied.

Yet another approach to the N-body problem is provided by special-purpose hardware like the GRAPE board (Ebisuzaki et al., 1993, Fukushige et al., 1991, Fukushige et al., 1996, Ito et al., 1991, Kawai et al., 2000, Makino, 1990, Makino and Funato, 1993, Makino et al., 1997, Okumura et al., 1993). It consists of custom chips that compute gravitational forces by the direct summation technique. By means of their enormous computational speed they can considerably extend the range where direct summation remains competitive with pure software solutions. A recent overview of the family of GRAPE-boards is given by Hut and Makino (1999). The newest generation of GRAPE technology, the GRAPE-6, will achieve a peak performance of up to 100 TFlops (Makino, 2000), allowing direct simulations of dense stellar systems with particle numbers approaching 10⁶. Using sophisticated algorithms, GRAPE may also be combined with P³M (Brieu et al., 1995) or tree algorithms (Athanassoula et al., 1998, Fukushige et al., 1991, Makino, 1991a) to maintain its high computational speed even for much larger particle numbers.

In recent years, collisionless dynamics has also been coupled to gas dynamics, allowing a more direct link to observable quantities. Traditionally, hydrodynamical simulations have usually employed some kind of mesh to represent the dynamical quantities of the fluid. While a particular strength of these codes is their ability to accurately resolve shocks, the mesh also imposes restrictions on the geometry of the problem, and onto the dynamic range of spatial scales that can be simulated. New adaptive mesh refinement codes (Klein et al., 1998, Norman and Bryan, 1998) have been developed to provide a solution to this problem.

In cosmological applications, it is often sufficient to describe the gas by smoothed particle hydrodynamics (SPH), as invented by Lucy (1977) and Gingold and Monaghan (1977). The particle-based SPH is extremely flexible in its ability to adapt to any given geometry. Moreover, its Lagrangian nature allows a locally changing resolution that ‘automatically’ follows the local mass density. This convenient feature helps to save computing time by focusing the computational effort on those regions that have the largest gas concentrations. Furthermore, SPH ties naturally into the N-body approach for self-gravity, and can be easily implemented in three dimensions.

These advantages have led a number of authors to develop SPH codes for applications in cosmology. Among them are TREESPH (Hernquist and Katz, 1989, Katz et al., 1996), GRAPESPH (Steinmetz, 1996), HYDRA (Couchman et al., 1995, Pearce and Couchman, 1997), and codes by Carraro et al. (1998), Davé et al. (1997), Evrard (1988), Hultman and Källander (1997), Navarro and White (1993). See Kang et al. (1994) and Frenk et al. (1999) for a comparison of many of these cosmological hydrodynamic codes.

In this paper we describe our simulation code GADGET (GAlaxies with Dark matter and Gas intEracT), which can be used both for studies of isolated self-gravitating systems including gas, or for cosmological N-body/SPH simulations. We have developed two versions of this code, a serial workstation version, and a version for massively parallel supercomputers with distributed memory. The workstation code uses either a tree algorithm for the self-gravity, or the special-purpose hardware GRAPE, if available. The parallel version works with a tree only. Note that in principle several GRAPE boards, each connected to a separate host computer, can be combined to work as a large parallel machine, but this possibility is not implemented in the parallel code yet. While the serial code largely follows known algorithmic techniques, we employ a novel parallelization strategy in the parallel version.

A particular emphasis of our work has been on the use of a time integration scheme with individual and adaptive particle timesteps, and on the elimination of sources of overhead both in the serial and parallel code under conditions of large dynamic range in timestep. Such conditions occur in dissipative gas-dynamical simulations of galaxy formation, but also in high-resolution simulations of cold dark matter. The code allows the usage of different timestep criteria and cell-opening criteria, and it can be comfortably applied to a wide range of applications, including cosmological simulations (with or without periodic boundaries), simulations of isolated or interacting galaxies, and studies of the intergalactic medium.

We thus think that GADGET is a very flexible code that avoids obvious intrinsic restrictions for the dynamic range of the problems that can be addressed with it. In this methods-paper, we describe the algorithmic choices made in GADGET  which we release in its parallel and serial versions on the internet,¹ hoping that it will be useful for people working on cosmological simulations, and that it will stimulate code development efforts and further code-sharing in the community.

This paper is structured as follows. In Section 2, we give a brief summary of the implemented physics. In Section 3, we discuss the computation of the gravitational force both with a tree algorithm, and with GRAPE. We then describe our specific implementation of SPH in Section 4, and we discuss our time integration scheme in Section 5. The parallelization of the code is described in Section 6, and tests of the code are presented in Section 7. Finally, we summarize in Section 8.