GADGET: a code for collisionless and gasdynamical cosmological simulations
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, P3M, and AP3M 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 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 P3M 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 106. Using sophisticated algorithms, GRAPE may also be combined with P3M (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,1 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.
Section snippets
Collisionless dynamics and gravity
Dark matter and stars are modeled as self-gravitating collisionless fluids, i.e., they fulfill the collisionless Boltzmann equation (CBE)where the self-consistent potential Φ is the solution of Poisson’s equationand is the mass density in single-particle phase-space. It is very difficult to solve this coupled system of equations directly with finite difference methods. Instead, we will follow the common N-body approach, where
Tree algorithm
An alternative to Fourier techniques, or to direct summation, are the so-called tree methods. In these schemes, the particles are arranged in a hierarchy of groups. When the force on a particular particle is computed, the force exerted by distant groups is approximated by their lowest multipole moments. In this way, the computational cost for a complete force evaluation can be reduced to order (Appel, 1985). The forces become more accurate if the multipole expansion is carried out to
Smoothed particle hydrodynamics
SPH is a powerful Lagrangian technique to solve hydrodynamical problems with an ease that is unmatched by grid based fluid solvers (see Monaghan, 1992, for an excellent review). In particular, SPH is very well suited for three-dimensional astrophysical problems that do not crucially rely on accurately resolved shock fronts.
Unlike other numerical approaches for hydrodynamics, the SPH equations do not take a unique form. Instead, many formally different versions of them can be derived.
Time integration
As a time integrator, we use a variant of the leapfrog involving an explicit prediction step. The latter is introduced to accommodate individual particle timesteps in the N-body scheme, as explained later on.
We start by describing the integrator for a single particle. First, a particle position at the middle of the timestep Δt is predicted according toand an acceleration based on this position is computed, viz.Then the particle is advanced according
Parallelization
Massively parallel computer systems with distributed memory have become increasingly popular recently. They can be thought of as a collection of workstations, connected by a fast communication network. This architecture promises large scalability for reasonable cost. Current state-of-the art machines of this type include the Cray T3E and IBM SP/2. It is an interesting development that ‘Beowulf’-type systems based on commodity hardware have started to offer floating point performance comparable
Tests of timestep criteria
Orbit integration in the collisionless systems of cosmological simulations is much less challenging than in highly non-linear systems like star clusters. As a consequence, all the simple timestep criteria discussed in Section 5.1 are capable of producing accurate results for sufficiently small settings of their tolerance parameters. However, some of these criteria may require a substantially higher computational effort to achieve a desired level of integration accuracy than others, and the
Discussion
We have presented the numerical algorithms of our code GADGET, designed as a flexible tool to study a wide range of problems in cosmology. Typical applications of GADGET can include interacting and colliding galaxies, star formation and feedback in the interstellar medium, formation of clusters of galaxies, or the formation of large-scale structure in the universe.
In fact, GADGET has already been used successfully in all of these areas. Using our code, Springel and White (1999) have studied the
Acknowledgements
We are grateful to Barbara Lanzoni, Bepi Tormen, and Simone Marri for their patience in working with earlier versions of GADGET. We thank Lars Hernquist, Martin White, Charles Coldwell, Jasjeet Bagla and Matthias Steinmetz for many helpful discussions on various algorithmic and numerical issues. We also thank Romeel Davé for making some of his test particle configurations available to us. We are indebted to the Rechenzentrum of the Max-Planck-Society in Garching for providing excellent support
References (99)
J. Comp. Phys.
(1995)J. Comp. Phys.
(1990)- et al.
NewA
(1997) NewA
(1996)- et al.
J. Comp. Phys.
(1974) - et al.
J. Comp. Phys.
(1987) - et al.
J. Comp. Phys.
(1999) - et al.
NewA
(1998) - et al.
Comput. Phys. Rep.
(1989) - et al.
J. Comp. Phys.
(1983)
NewA
J. Comp. Phys.
Comp. Phys. Comm.
MNRAS
ApJ
SIAM
J. Sci. Stat. Comp.
MNRAS
MNRAS
The use of supercomputers in stellar dynamics
Nature
ApJS
Comput. Phys.
ApJS
ApJ
MNRAS
ApJ
ApJ
ApJ
PASJ
ApJS
MNRAS
ApJ
PASJ
ApJ
MNRAS
ApJS
ApJ
A&A
Computer Simulation Using Particles
ApJ
ApJ
A&A
Science
ApJ
PASJ
Cited by (1701)
Fast forward modelling of galaxy spatial and statistical distributions
2024, Journal of Cosmology and Astroparticle PhysicsThe long-lasting effect of X-ray pre-heating in the post-reionization intergalactic medium
2024, Monthly Notices of the Royal Astronomical Society