We provide introductory explanations and illustrations of the N-body hydrodynamics code Charm N-body GrAvity solver (ChaNGa). ChaNGa simulates the gravitational motion and gas dynamics of matter in space with the goal of modeling galactic and/or cosmological structure and evolution. We discuss the algorithm for leapfrog integration and smoothed particle hydrodynamics and computer science concepts used by the program, including the binary data structure for the particle positions. Our presentation borrows from the doctoral dissertation of J. Stadel, U. Washington, 2001. Problems are provided in order to use ChaNGa to learn or solidify some cosmological concepts.
I. INTRODUCTION
N-body hydrodynamic simulations of the cosmos and isolated galaxies1,2 have been crucial in establishing Lambda cold dark matter (Lambda refers to dark energy) as the standard model of cosmology.3 Due to the fundamental nature of gravity, computationally challenging densities occur in astrophysical situations. Simulations developed by several groups have yielded trustworthy predictions of the incredibly diverse structure of the cosmos and galaxies, and several simulations of very differing kinds have given the same results, thanks to the AGORA project.4
ChaNGa (Charm N-body GrAvity solver) is an example of a highly parallel N-body hydrodynamic code.5 The physics of star formation, supernovae feedback, and cooling are included in ChaNGa.
Our goal is to provide an accessible, but thorough explanation of how to use this N-body hydrodynamic code by adding to ChaNGa's GitHub wiki6 and many YouTube presentations. We also hope to provide a deeper appreciation for how computer science data structures are used to construct ChaNGa. For example, ChaNGa is written in the parallel programming language CHARM++, which won a Gordon Bell prize.7 Although we have not attempted to provide an exhaustive description of the code, we focus on the details that are most relevant to the physics.
Sections II–IV discuss gravitational motion and the ChaNGa's approach to treat it. Section V describes several considerations for a simulation of an infinite (or at least very large) universe. Section VI describes a technique for softening gravitational forces at small distances to avoid problems that would otherwise arise in the calculation of the forces. Section VII describes the method by which ChaNGa accounts for the presence of gas in the universe. Section VIII gives some of the results of ChaNGa, and Sec. IX provides some suggested problems.
II. NEWTONIAN GRAVITY AND NUMERICAL INTEGRATION
Gravity is responsible for Kepler's three laws, which are commonly derived in junior-level mechanics without any mention of computers. The difference between the N = 2 and N > 2 systems is huge despite the fact that Newton's laws and the universal law of gravitation are the same for both.
There are several reasons why the kick-drift-kick form of the leapfrog algorithm is favored by ChaNGa. Paramount among these is a concept known as multi-stepping.5 The number of computations in Eqs. (5)–(7) made during a fixed time interval largely determines the accuracy of this integration algorithm. Because the method is an approximation in which particles move with a fixed velocity between time steps, the smaller the time step, the more accurate the trajectories. There is a tradeoff in terms of run-time, and so rather than pursuing increasingly small time steps, the goal is usually to find the number of time steps compared to the duration of the time interval for which the accuracy is good enough. Because gravity is spatially dependent, what constitutes sufficient temporal resolution is closely linked with the particle density. Regions with higher particle density will naturally require smaller time steps to achieve the same level of physical accuracy that can be attained in a lower density region with a larger time step.1,9,11
One approach to this problem would be to increase the number of time steps globally to be suitable to the desired level of accuracy in high density regions, but doing so would introduce a large number of extraneous calculations for low density regions and, thus, would be computationally wasteful.1 The solution is to assign time steps on a per-particle basis; that is, each particle gets its own time step in accordance with the density of its surrounding particles (multi-stepping). The kick-drift-kick method ensures that such an individualized time-step assignment is possible, as long as the time step sizes of all particles are related by powers of two.9,12
III. A HIERARCHICAL DATA STRUCTURE FOR THE PARTICLE POSITIONS: THE BARNES–HUT ALGORITHM
Due to the complex nature of the N-body problem, another challenge is run-time efficiency. To understand the challenge, note that the net acceleration in the leapfrog algorithm in Eqs. (5)–(7) is the result of the large sum of gravitational interactions in Eq. (2) for . From the form of Eq. (2), we see that for a system of N particles, each particle experiences N – 1 gravitational interactions. If we were to write a program with just the information from Sec. II, our program would need to perform N – 1 force calculations N times or calculations at each time step.8 This number might be acceptable for small-scale computations, say for N = 100–1000, but the current state of computational astrophysics is such that simulations involve N well in the millions and often evolve over thousands of time steps.9,11 Programs need to accommodate this level of scale, which requires methods to speed up the gravitational force calculation process.
The solution employed by ChaNGa is to use a data structure designed to speed up the access to the huge amount of initial positions data: A tree structure. The most well-known implementation of such a structure is the Barnes–Hut algorithm.13 Barnes–Hut offers an approximation of the gravitational force acting on a particle in logarithmic time ( ], making the total force calculation per time step , a gargantuan improvement from the polynomial time direct summation method.
In the following, we introduce the standard Barnes–Hut algorithm to familiarize the reader with the general approach to the Gordon-Bell-prize-winning solution for speeding up the N-body problem. We will then discuss some of the modifications employed by ChaNGa and the motivations behind them. Readers who wish a more in-depth or specific treatment of Barnes–Hut might be interested in the seminal 1986 paper,13 or more accessible sources such as Refs. 15 and 9.
The fundamental principle behind the Barnes–Hut algorithm is to reduce the number of force calculations performed for each particle by applying a center-of-mass approximation for the force exerted by collections or clusters of distant particles.17 More explicitly, for distant clusters—particles that are far from the particle whose acceleration we want to calculate but are close to each other17—the algorithm determines the center of mass and total mass of the cluster and then performs a single gravitational force calculation instead of many force calculations, which contribute little to improve the accuracy of the total force.
As noted in Ref. 9, the nature of gravitation is such that for a sufficiently distant cluster of identical particles, the force contribution associated with each particle will be roughly the same (both in magnitude and direction), so the approximation is justified. Of course, the distance constitutes the main underlying assumption for the approximation, meaning it does not hold for nearby particles, for which direct summation is still performed. Because this approach relies heavily on distance, we need to establish a way to determine what constitutes “far enough” for a center-of-mass approximation and “close enough” for particles to be considered part of the same cluster. We might be tempted to simply compare the distances of all the particles in the simulation, but such a process would be ,9 and be no better than direct summation.
The Barnes–Hut algorithm is a clever solution to this problem, which starts with a process known as domain decomposition, or the repeated decomposition of the initial volume into smaller subvolumes. It is useful to define the widely used term node in the context of hierarchical data structures: “each node represents a portion of the 3D space containing the particles in that volume.”5 The Barnes–Hut algorithm uses an octree decomposition data structure,9,13 which involves recursively dividing cubical volumes into eight subvolumes, hence, the name. (This decomposition is not the scheme ChaNGa uses, but understanding the basic octree decomposition method is helpful as an introduction to the method in ChaNGa.) A summary of octree decomposition is as follows:
-
Begin with the volume containing all particles, this volume is the root node of the octree.
-
Make the first decomopostion by dividing the root node evenly into eight subvolumes.
-
If there are no particles in one of the subvolumes/nodes, discard it.
-
If two or more particles are present in a subvolume, divide it evenly into eight daughter subvolumes/nodes. Then inspect each daughter node individually.
-
Repeat steps 2–4 until each node contains at most one particle.
In total, the complexity of constructing the octree is .13
Figure 1 gives an example of this process for a system of four particles. The nodes vary in size according to their corresponding level in the tree structure. Only boxes containing more than one particle are subdivided into daughter nodes, and division is performed recursively until each node contains one or zero particles. In practice, the number of particles which is the threshold for further decomposition is higher than one and is a parameter that can be varied to adjust the resolution. See Ref. 17 for an excellent description of this process.
A visual representation of an Octree decomposition for four particles taken from Ref. 14. The figure on underneath gives an alternative visualization of the decomposition. The purple node at the top corresponds to the entire original volume. The next layer down is orange-pink with three kinds of boxes in that row: no color for no particles, orange–pink with a black dot indicating a single particle, and an orange–pink box with no black dots indicating that there are two or more particles in that node (two particles in this example). These nodes require further decomposition.
A visual representation of an Octree decomposition for four particles taken from Ref. 14. The figure on underneath gives an alternative visualization of the decomposition. The purple node at the top corresponds to the entire original volume. The next layer down is orange-pink with three kinds of boxes in that row: no color for no particles, orange–pink with a black dot indicating a single particle, and an orange–pink box with no black dots indicating that there are two or more particles in that node (two particles in this example). These nodes require further decomposition.
We still need to discuss how this hierarchical data structure expedites the force calculation, an insight which may yet elude all but the most computer-savvy readers.
The next step involves the center of mass calculations. The implementation of the algorithm as described in Ref. 13 is such that the data structure representing a node in the tree has attributes corresponding to the total mass and center of mass for the particles which it contains. For nodes containing only one particle (“leaf nodes” in the language of the tree metaphor), these values are trivial, so it is most efficient to propagate this information backward through the data structure, i.e., leaves-to-root, a process which is also .13
If we examine the analogous two-dimensional Quadtree data structure depicted in Fig. 2, we see that the center of mass and total mass calculation would begin with node 4, whose attributes would be calculated using the positions and masses of particles a and b. We can then use these attributes of node 4, as well as the mass and position of particle d, to compute the total mass and center of mass of node 2 and so on. Nodes containing more than one particle can then be treated as pseudoparticles with their own associated mass (the total mass of the particles) and position (the center of mass position).
A two-dimensional analogue of the Barnes–Hut Octree, known as a Quadtree. Such figures are better understood initially by removing all the subvolumes and starting from scratch. Only domain decomposed subvolumes that have more than one particle are shown in (a). The numbers next to the circles in (b) refer to a particular node in the text. This image was taken from Ref. 16.
A two-dimensional analogue of the Barnes–Hut Octree, known as a Quadtree. Such figures are better understood initially by removing all the subvolumes and starting from scratch. Only domain decomposed subvolumes that have more than one particle are shown in (a). The numbers next to the circles in (b) refer to a particular node in the text. This image was taken from Ref. 16.
Once all the pseudoparticle nodes in the system have been assigned their mass and center of mass attributes, all the machinery necessary to begin force calculations is in place. The algorithm for calculating the force on a given particle p involves a traversal of the Octree, where is the side length of the cubical region represented by the current node, D is the distance between p and the node's center of mass, and θ is an accuracy parameter set at the start of the simulation, usually . The algorithm for each particle is13
-
Start at the root node.
-
If , compute the gravitational force between the current node and p, and add it to the total force acting on p.
-
Otherwise, traverse one layer down the tree (away from the root node) and perform step 2 for each daughter node.
Once these steps have been completed, the gravitational influence of each particle on p is accounted for either by a direct force calculation or by a center of mass approximation. For large N, this process involves performing order force calculations for each particle, and, thus, the overall run time of the algorithm is .
IV. BARNES–HUT MODIFICATIONS
The description we have provided of the Barnes–Hut gravitational force computation closely follows Ref. 13, but as mentioned, there are several changes in the ChaNGa code. Most of these changes are described in Ref. 20. One notable change is that ChaNGa does not actually use the center of mass of faraway clusters to approximate their gravitational influence. Instead, ChaNGa performs an operation known as a multipole expansion (usually discussed for electric fields) for improved force accuracy,5 see Ref. 17 for an excellent description of multipole expansions for N-body simulations. Also see the discussion of multipole expansions in Ref. 18. Due to the similarities between Newtonian gravity and the electric field, multipole expansions can be used to approximate the gravitational potential produced by clusters of massive particles.
A multipole expansion effectively amounts to an infinite sum of terms with increasing angular dependence.19 The first (or zeroth order) term is a monopole, or a point source, which has no angular dependence. The first few higher order terms are the dipole, quadrupole, and hexadecapole terms in the electromagnetism context,18,19 or in Ref. 20 in the N-body simulation context. ChaNGa performs a multipole expansion to third order, including all terms up to and including hexadecapole,5 which is faster and more accurate than a quadrupole order expansion.20
Another notable difference between the Barnes–Hut algorithm and that employed by ChaNGa is the precise structure of the spatial decomposition tree. As we discussed in Sec. III, the algorithm recursively divides the spatial domain into even sections of eight, resulting in a data structure called an octree (see Fig. 1). Instead, the spatial decomposition performed by ChaNGa employs a binary tree rather than an octree, meaning the domain is recursively divided into sections of two rather than eight. This change is consistent with the approach used by the N-body code PKDGRAV,20 from which ChaNGa inherits much of its gravitational force calculation.5
The justification for using a binary tree rather than an octree is that a binary tree structure offers advantages for parallelization, particularly when distributing the computation over an arbitrary number of processors,20 a feature which has become of central importance in high performance computing.
Although the precise details of ChaNGa's spatial binary tree implementation are not easy to learn from the literature, ChaNGa explicitly inherits much of its gravitational force calculation from PKDGRAV,5 and the details of which can be found in Ref. 20. Because of this intellectual inheritance, a brief discussion on PKDGRAV's domain decomposition method is warranted to provide a better understanding of what is implemented in ChaNGa. Aside from the overall difference in an organizational structure, a crucial difference between the Barnes–Hut Octree and the binary trees used by contemporary simulation codes is that in binary tree data structures space is often not divided up evenly. Instead, daughter nodes are sized dynamically in accordance with some bisection scheme, often related to the number of particles per daughter node. An example of such a binary tree algorithm (in two dimensions) is the k-D tree, depicted in Fig. 3(a). This algorithm involves recursively bisecting nodes through their longest axis such that both daughter nodes contain approximately the same number of particles. Note that in Fig. 3, the first bisection is done such that each daughter node contains 15 particles, then 8, and finally 4.
Images representing a standard k-D tree domain decomposition (a) and the spatial binary tree decomposition used by PKDGRAV (b). The k-D tree is constructed with a maximum of four particles per leaf cell. On the right, the bounds have been squeezed before bisection. The dashed lines indicate bisectors. The differences between the two decompositions are much more pronounced for a more highly clustered distribution. In that case, the spatial locality of the spatial binary tree would be much improved over the k-D tree. k-D trees can have pathological problems with certain distributions (Ref. 20), so the spatial binary tree is the best choice. The name spatial binary tree comes from “spatially bisecting the bounding box.” This image and caption are taken from Ref. 20.
Images representing a standard k-D tree domain decomposition (a) and the spatial binary tree decomposition used by PKDGRAV (b). The k-D tree is constructed with a maximum of four particles per leaf cell. On the right, the bounds have been squeezed before bisection. The dashed lines indicate bisectors. The differences between the two decompositions are much more pronounced for a more highly clustered distribution. In that case, the spatial locality of the spatial binary tree would be much improved over the k-D tree. k-D trees can have pathological problems with certain distributions (Ref. 20), so the spatial binary tree is the best choice. The name spatial binary tree comes from “spatially bisecting the bounding box.” This image and caption are taken from Ref. 20.
The k-D tree offers a good first look into spatial binary trees with dynamically sized nodes, but it also presents problems in terms of force error and run time efficiency.20 The decomposition method used by PKDGRAV is depicted in two dimensions in Fig. 3(b). By comparing this diagram to its neighbor representing the k-D tree, it is clear that there is more going on here than explained in the caption (Fig. 4). You can gain a better understanding of spatial binary trees by reading the appropriate section of Ref. 20.
The left-hand side is a square divided into four sub-squares with the “curve” going through each of the four sub-squares. The curve looks like the letter “N,” which might not look like a curve. The figure is helpful in understanding the meaning of a space-filling curve, because we can see the interpretation of the “space” that is required to understand the acronym “space-filling curve.” That is, there are only four “points”—instead of an infinite number—in this space. The right-hand figure is intended to help generalize the point of the left-hand subfigure. The figures are from Ref. 22.
The left-hand side is a square divided into four sub-squares with the “curve” going through each of the four sub-squares. The curve looks like the letter “N,” which might not look like a curve. The figure is helpful in understanding the meaning of a space-filling curve, because we can see the interpretation of the “space” that is required to understand the acronym “space-filling curve.” That is, there are only four “points”—instead of an infinite number—in this space. The right-hand figure is intended to help generalize the point of the left-hand subfigure. The figures are from Ref. 22.
To better understand what is shown in Fig. 3, we note that the spatial binary tree decomposition really has only one fundamental difference from the k-D tree: The “squeezing” of daughter nodes to minimize the volume they represent. At every step, PKDGRAV's binary tree compresses the volume of the node being examined to the smallest rectangle (rectangular prism in three dimensions) that contains all the particles in the node. Because the bounding box is now considered to be the node's volume, it is bisected through its longest axis into two daughter nodes containing roughly equal numbers of particles, which are then squeezed themselves. This process is repeated until each node is under the maximum particle count for a leaf node, which for ChaNGa is usually around 8–12.5
We have so far discussed two of the fundamental modifications to the Barnes–Hut algorithm employed by ChaNGa, namely, the use of a hexadecapole order multipole expansion for the gravitational force approximation at large distances and the use of a spatial binary tree structure rather than an octree. There is another aspect of the gravitational force calculation that should be mentioned, namely, ChaNGa's approach to parallelization, because large-scale parallelization is a primary reason for the code.21 A detailed description of the parallelization process can be found in Sec. 4 of Ref. 5.
We have described how ChaNGa uses a spatial binary tree to handle the gravitational force calculation. Although it is correct that the binary tree is the fundamental data structure, ChaNGa's force calculation involves not one but a number of spatial binary trees, each representing a subregion of the simulation volume containing a subset of the overall particle count, divided among the processors allotted to the computation.5 To facilitate this approach, ChaNGa performs an initial domain decomposition step not described by the spatial binary tree or octree, in which the simulation volume is divided into a number of subregions containing an equal number of particles according to a space-filling curve algorithm.5 The space-filling curve passes through each discrete spatial cell containing a particle, Fig. 4. The resulting spatial decomposition looks different from a binary tree or octree. This space-filling curve algorithm initially occurs to balance the computational load across the processors being used. The code iterates through a number of potential decompositions until all bins (subregions) are sufficiently optimal.5 All the particles in each bin are then assigned to a data object called a tree piece, so that each tree piece represents a subset of the overall volume. Each tree piece then performs a spatial binary decomposition of its assigned subvolume (using a bounding box method similar or identical to the one we have described), and the gravitational force calculation begins.5
The term “tree piece” is closely associated with another concept in ChaNGa; a chare. Chares are important in forming checkpoints for long simulations, and chares are essentially tree pieces.5 Note that, because particles are assigned to tree pieces in accordance with the initial space-filling curve decomposition, no single processor has direct access to all the particles in the simulation, but all particles interact with each other gravitationally. Because of this dichotomy, ChaNGa also facilitates communication across processors to allow remote access to nodes that are part of different tree pieces.5 Communication across processors involves a network.
V. BOUNDARY CONDITIONS AND EWALD SUMS
We know that the observable universe spans a distance of approximately 13 × 109 light years. Even very large simulations, such as those performed as part of the IllustrisTNG project,23 have computational volumes with widths only in hundreds of millions of light years. Thus, even exceptionally large simulations model regions, which constitute only a small fraction of the universe. Because gravity acts at large distances, this limitation is a problem for realistic simulations. Before discussing an approach to this problem, we first discuss an important concept of the currently accepted cosmological model: The cosmological principle. The cosmological principle states that the universe is both isotropic and homogeneous.24 Being isotropic means that the universe looks the same in all directions. Similarly, homogeneous means that the universe looks the same at all locations. Readers who have looked at the night sky may object, because the universe does not look the same in all directions and there are distinct features such as stars and constellations. Similarly, others might argue that the universe must look a little different for someone in the Andromeda galaxy. The cosmological principle needs to be understood in the context of distance scales: On some scales, the universe is isotropic and homogeneous, but on others, it definitely is not. The cosmological principle does not state that on a scale such as our observational perspective, the universe cannot contain distinct features, but rather that at sufficiently large scales (relative to our perspective), the two assumptions hold.25
The cosmological principle provides justification for modeling a simulation volume as a single cell in an infinite (or at least very large) three-dimensional grid of perfectly identical cubes, which is the approach taken by ChaNGa. Such codes use periodic boundary conditions, because the structure of the simulation volume is repeated over a periodic lattice.20 This approach addresses the issues that arise from a finite simulation volume, but it is then necessary to treat the gravitational force calculation over an infinite lattice, which might seem no less daunting. The solution adopted by ChaNGa is to break up the calculation into long-range and short-range components.11,20 The short-range calculation is performed as an extension of the Barnes–Hut algorithm described in Sec. III by including a number of neighboring lattice cells (usually 26 of them20) in the Barnes–Hut force calculation. The long-range gravitational contribution is then accounted by Ewald sums,5,20,27 where the simulation cell is repeated exactly in every direction so that the long range force can be calculated as a sum over all these repeated cells.
VI. FORCE SOFTENING
Another issue related to the gravitational force calculation is how to handle the gravitational force between particles at very small distances. From Eq. (1), we can see that the interparticle force increases very quickly at small particle separations. This infinite force for zero separation presents a problem, because it can be difficult to handle very large forces computationally,9 leading to unphysical results.20 Hence, it is important for an N-body code to implement some kind of softening to impose a limit on the magnitude of the force. In ChaNGa, this force softening is closely related to the handling of smoothed particle hydrodynamics (see Sec. VII) and involves a spline softening kernel,5,20 which effectively cuts off the gravitational force at zero interparticle separation, while maintaining the usual Newtonian gravity at distances greater than the softening length.26
VII. SMOOTH PARTICLE HYDRODYNAMICS
We see that a robust and efficient computational treatment of gravity is no simple task. There remain several nuances, which have been discussed only superficially or have been omitted entirely.5,20,26,28 We could write a book just on the workings of a single sufficiently advanced N-body code.20 Nevertheless, there is another crucial element of ChaNGa that we would be remiss to not discuss at least briefly: gas dynamics.
Everything we have discussed about ChaNGa so far has been related to predicting the motion of infinitesimal point masses under the influence of a mutual gravitational force. Clearly, this calculation is crucial for computational astrophysics: Gravity is a dominant force in the cosmos. Nevertheless, simulating the motion of particles under the influence of gravity would be sufficient only if we lived in a universe that consisted purely of discrete, massive bodies such as stars and dark matter. In particular, a purely gravitational simulation would ignore the significant amounts of gas and dust present in our universe, which also play a significant role in galaxy formation and cosmological structure.9 Gravity influences the motion of gas and dust, but additional considerations need to be made for the dynamics of the gas itself. To handle gas physics, N-body codes turn to fluid mechanics and treat the gas as a continuous fluid medium, a reasonable approximation at a macroscopic scale.9,29 The fundamental equation of fluid dynamics is the nonlinear Navier–Stokes equation. This nonlinearity leads to many more challenges.
There exist two approaches to fluid dynamics: the Eulerian description and the Lagrangian description.29 The Eulerian description considers what is happening at fixed locations (or points) in space as a fluid medium flows through these points. In contrast, the Lagrangian description considers the dynamic evolution of individual pieces of matter representing “fluid elements” as they travel through space.29 Broadly speaking, there are two types of N-body simulation codes that treat gas dynamics: grid-based (also known as mesh-based) codes, which follow the Eulerian description by considering a simulation grid and tracking the flow of gas through that grid, and particle-based codes, such as ChaNGa, which follow the Lagrangian description by tracking the spatial motion of individualized gas “parcels.”26
One of the advantages to the Lagrangian approach is that it is a fairly natural extension of the Barnes–Hut based gravitational force calculation described in Sec. III.5 The particles involved in the gravitational force calculation are also used to discretize a gas, effectively representing the “parcels” for the Lagrangian approach.9 This important and elegant step is accomplished through a process known as smoothed particle hydrodynamics, which uses the particles present in the simulation to derive continuous fluid quantities (such as the pressure and temperature) spanning the surrounding region of space.9,11 In this sense, the particles can be thought of as being “smoothed out” over space, smudged into a continuous particle-fluid that represents the gas.
(Color online) Depiction of a smoothing kernel W applied to a discrete set of particles. The blue-filled circles represent the particles, h is the smoothing length, and the curve represents the value of the smoothing kernel for all points in space. The frame of reference in the lower left is intended to help understand the arguments of W. The vector is shown in the plane of the blue-filled circles (Ref. 30).
(Color online) Depiction of a smoothing kernel W applied to a discrete set of particles. The blue-filled circles represent the particles, h is the smoothing length, and the curve represents the value of the smoothing kernel for all points in space. The frame of reference in the lower left is intended to help understand the arguments of W. The vector is shown in the plane of the blue-filled circles (Ref. 30).
VIII. RESULTS
So far we have introduced N-body hydrodynamic simulation methods at the conceptual level, because most researchers do not write code from scratch. In the following, we give a brief overview of the user's perspective when running a simulation in ChaNGa on a laptop or desktop computer. Running a ChaNGa simulation on a supercomputer requires considerably more infrastructure.31
Given that you have access to a computer running Linux, doing a simulation in ChaNGa requires only the code itself, a parameter file, and an initial conditions file. ChaNGa is a complicated piece of software, and as such is supported by a number of dependencies. This complexity, as well as that of its dependencies, can lead to some difficulties in the installation process. For those having trouble installing the software, see Ref. 32 and the supplementary material, which provides some detailed notes on ChaNGa installation.33
ChaNGa simulates the gravitational motion and gas dynamics of matter in space with the goal of modeling galactic and/or cosmological structures. The initial conditions' file is a binary file (not human readable) that provides the initial positions and velocities of all particles in the simulation.11 These initial conditions cannot be arbitrarily assigned and must be computed in a manner that approximates observations of the early universe via the cosmic microwave background measured by the Wilkinson microwave anisotropy and Planck probes.34,35 For a discussion of the generation of initial conditions for the AGORA project using MUSIC,36 see Chap. 1 of Ref. 11. A more elementary and well-documented program for producing initial conditions that works well with ChaNGa is pyICs.37 The initial conditions provide the starting point by giving the initial position and velocity values for all the particles; these values are assigned to be approximately consistent with observation.
We used the software visualization tool pynbody40 to generate the image in Fig. 6 based on the cube300 simulation.42 Pynbody is suggested because it is particularly ChaNGa friendly, and the documentation is excellent. The fibrils of density corresponding to light blue are of great importance for transporting gas to and from galaxies, and their physics is of great interest.41 The purple dots in the filament represent galaxies.
A visualization of the cosmos at redshift zero based on a cube300 test simulation in the ChaNGa directory. (Reference 42 provides details for running cube300.) Note the web-like pattern of the density: In many regions, the density is zero. The fibril-like paths at densities indicated by light blue are often called filaments.
A visualization of the cosmos at redshift zero based on a cube300 test simulation in the ChaNGa directory. (Reference 42 provides details for running cube300.) Note the web-like pattern of the density: In many regions, the density is zero. The fibril-like paths at densities indicated by light blue are often called filaments.
Our goal has been to extend the presentations of the N-body hydrodynamic code ChaNGa given in Refs. 5 and 31 to make using ChaNGa more accessible. Some advanced topics at the heart of ChaNGa have regretfully been left out. It is hoped that others with experience using ChaNGa might take up the challenge of providing a part two on possible topics such as the CHARM++ language and how it works with various kinds of hardware.
IX. SUGGESTED PROBLEMS
A few suggested problems are provided to guide the reader in both the analysis of the simulation results and basic concepts in cosmological physics.
Problem 1. Units in ChaNGa are interesting and require attention. Provide a detailed derivation of the results for mass units in ChaNGa wiki example number two under “Understanding Units,” which is under “Code Units.”38 The derivation will get the reader ready to understand the units and value of the Hubble constant in problem2. Reference 39 is extensive and requires some study before finding “Code Units.”
Problem 2. Imagine that you are responsible for creation and that one of your tasks is to create a multiverse. Simulate a universe that has a Hubble constant H0 that is, say, roughly half the value of the presently accepted value. Use ChaNGa and cube300 to generate a universe and demonstrate that the simulation behaves as expected. Note: This problem is a bit—possibly very—tricky, but once the data are found/obtained and compared, the result is very clear.
This problem depends strongly on understanding the material in the supplementary material.42 For example, it provides the Hubble constant parameter and explains the issue of initial conditions for this simulation.
Problem 3. Use Pynbody to demonstrate that, in fact, cube300 is a dark matter only simulation. The installation documentation for Pynbody from GitHub is quite good.41 Some knowledge of Python 3 is necessary to use Jupyter to access the low-level abilities of Pynbody for the solution to this problem.
ACKNOWLEDGMENTS
The authors gratefully acknowledge Reed College for summer research internship support. J.W.P. would also like to acknowledge ACCESS Allocation (No. AST200020) for supporting this project.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts of interest to disclose.