9. Comparison with GADGET-2

In this section we shall perform equivalent simulations with CONCEPT and the well-known GADGET-2 code. Due to the two codes utilising different numerical schemes, as well as the fact that CONCEPT — unlike GADGET-2 — is designed to be consistent with general relativistic perturbation theory by default, some subtleties need to be taken into account in order to perform simulations with CONCEPT that are truly equivalent to those of GADGET-2.

If you are unfamiliar with GADGET-2 or have no interest in running GADGET-2 equivalent simulations with CONCEPT, you may skip this section. However, many of the tricks to come are not exclusively related to GADGET-2.

CON CEPT

As always, we start with a CONCEPT parameter file, which you should save as e.g. param/tutorial-9:

param/tutorial-9

# Input/output
if _gen:
    initial_conditions = {        'name'   : 'GADGET halo',        'species': 'matter',
        'N'      : _size**3,
    }
else:
    initial_conditions = f'{path.ic_dir}/{param}_a={a_begin}'
output_dirs = {
    'snapshot': path.ic_dir if _gen else f'{path.output_dir}/{param}',
    'powerspec': ...,
}
output_bases = {
   'snapshot': param if _gen else 'snapshot',
}
output_times = {
    'snapshot' : a_begin if _gen else '',
    'powerspec': '' if _gen else [0.1, 0.3, 1.0],
}snapshot_type = 'gadget'
# Numerics
boxsize = 512*Mpc/h
potential_options = 2*_size

# Cosmology
H0 = 67*km/(s*Mpc)
Ωb = 0.049
Ωcdm = 0.27
a_begin = 0.02

# Physicsrealization_options = {    'gauge'     : 'synchronous',    'back-scale': True,}# Non-parameter helper variables
_size = 64_gen = False

To perform the same simulation with both CONCEPT and GADGET-2, we shall want to start from a common initial snapshot. The parameter file is set up to produce this initial snapshot when _gen is True. To create the snapshot then, do

./concept \
    -p param/tutorial-9 \
    -c "_gen = True"

(as always, feel free to tack on the number of processes you want with -n). The snapshot will be dumped to the ic directory.

As we have snapshot_type = 'gadget', the snapshot will be in GADGET format. Furthermore, naming the matter component 'GADGET halo' ensures that this component gets mapped to GADGET particle type 1 (conventionally used for cold dark matter, see table 3 in the user guide for GADGET-2 ).

Starting from this snapshot we can now perform the CONCEPT simulation by running with -c "_gen = False". As this is also the value set in the parameter file, we may in fact start the simulation simply by

./concept -p param/tutorial-9

Note that the initial_conditions is set to the path of the generated snapshot when _gen is False.

The power spectrum outputs of the CONCEPT simulation will be dumped to output/tutorial-9 as usual.

GADGET-2

For GADGET-2, the first thing we need is the source code itself. If you have installed CONCEPT the easy way, you already have it! Let’s make a copy just for use with this tutorial:

(source concept && cp -r "${Gadget2_dir}" "${output_dir}"/tutorial-9/)

You now have the complete GADGET-2 code in output/tutorial-9/Gadget2.

The Makefile of GADGET-2 needs to be set up with correct path information for its dependencies. Furthermore, various options need to be set in order for the GADGET-2 simulation to come to be equivalent to the CONCEPT simulation. Last but not least, we need a GADGET-2 parameter file equivalent to the CONCEPT parameter file. All of this can be conveniently achieved using the gadget utility included with CONCEPT:

./concept \
    -u gadget output/tutorial-9/Gadget2 \
    -p param/tutorial-9

The output/tutorial-9/Gadget2 directory now has a properly set up Makefile and a parameter file called param. The output times of the original parameter file have also been copied to an outputlist file, similarly placed in output/tutorial-9/Gadget2.

Note

The parameters specified in the GADGET-2 Makefile and param file include the cosmology, the grid size of the PM grid, the particle softening lengths and more, all set to match those used by the specific CONCEPT simulation as specified by the parameter file together with default CONCEPT values.

The GADGET-2 TreePM gravitational method is used, which is quite similar in essence to the P³M method used by the CONCEPT simulation, the difference being that the tree approximates the gravitational short-range force by grouping particles together.

Lastly, the values of physical constants employed by CONCEPT and GADGET-2 differ slightly, which will make a small difference for the comparison. If the GADGET-2 source code you are using is the one that goes with the CONCEPT installation, all such constants of the GADGET-2 source code have been patched to match the values used by CONCEPT.

With the various GADGET-2 files in place, let’s build the code and launch the simulation:

(source concept \
    && cd output/tutorial-9/Gadget2 \
    && make clean \
    && make \
    && $mpiexec -n 4 ./Gadget2 param)

where you may wish to change the specified number of processes deployed.

As GADGET-2 is not able to produce power spectra, the results from the simulation are raw snapshots, placed along with other files in output/tutorial-9/Gadget2/output. To investigate these snapshots, we can make use of the info utility:

./concept -u info output/tutorial-9/Gadget2/output

You should find that the snapshots indeed have the same values of the scale factor \(a\) as specified in the parameter file (perhaps with small numerical errors).

To measure the power spectra of the particle distributions contained within the snapshots, we may make use of the powerspec utility:

./concept \
    -u powerspec output/tutorial-9/Gadget2/output \
    -p param/tutorial-9

(as always, you may choose to throw more processes at the task with -n).

Comparison and adjustments

You should now have CONCEPT power spectra in output/tutorial-9 and corresponding GADGET-2 power spectra in /output/tutorial-9/Gadget2/output. We can of course look at each pair of plots, but for a proper comparison we should plot the power spectra together in a single plot. You may use the script below:

output/tutorial-9/plot.py

import glob, os, re
import numpy as np
import matplotlib; matplotlib.use('agg')
import matplotlib.pyplot as plt

# Particle Nyquist frequency
Mpc = 1
h = 0.67
boxsize = 512*Mpc/h
num_particles = 64**3
k_nyquist = 2*np.pi/boxsize*np.cbrt(num_particles)/2

# Read in CO𝘕CEPT data
def read(dirname):
    P = {}
    for filename in sorted(glob.glob(f'{dirname}/powerspec*'), key=os.path.getmtime):
        if filename.endswith('.png'):
            continue
        with open(filename, mode='r', encoding='utf-8') as f:
            header = f.readline()
        a = float(re.search('a = (.+)', header).group(1).rstrip('.'))
        k, P[float(f'{a:.3f}')] = np.loadtxt(filename, usecols=(0, 2), unpack=True)
    return k, P
this_dir = os.path.dirname(os.path.realpath(__file__))
k, P_concept = read(this_dir)

# Read in GADGET-2 data
k, P_gadget = read(f'{this_dir}/Gadget2/output')

# Plot
fig, axes = plt.subplots(2, sharex=True)
for a, P in P_concept.items():
    axes[0].loglog(k, P, label=f'$a = {a}$ (CO$N$CEPT)')
    if a in P_gadget:
        axes[0].loglog(k, P_gadget[a], 'k--')
        axes[1].semilogx(k, 100*(P/P_gadget[a] - 1))
for ax in axes:
    ylim = ax.get_ylim()
    ax.plot([k_nyquist]*2, ylim, ':', color='grey')
    ax.set_ylim(ylim)
ax.text(
    k_nyquist,
    ylim[0] + 0.15*np.diff(ylim),
    r'$k_{\mathrm{Nyquist}}\rightarrow$',
    ha='right',
)
axes[0].set_xlim(k[0], k[-1])
axes[1].set_xlabel(r'$k\, [\mathrm{Mpc}^{-1}]$')
axes[0].set_ylabel(r'$P\, [\mathrm{Mpc}^3]$')
axes[1].set_ylabel(r'$P_{\mathrm{CO}N\mathrm{CEPT}}/P_{\mathrm{GADGET}} - 1\, [\%]$')
axes[0].legend()
axes[0].tick_params('x', direction='inout', which='both')
axes[1].set_zorder(-1)
fig.tight_layout()
fig.subplots_adjust(hspace=0)
fig.savefig(f'{this_dir}/plot.png', dpi=150)

Store the script as e.g. output/tutorial-9/plot.py and run it using

./concept -m output/tutorial-9/plot.py

The upper subplot — with absolute power spectra — of the generated output/tutorial-9/plot.png should show a good qualitative match between the two codes, whereas the lower subplot — with relative power spectra — should reveal a disagreement of several percent. This disagreement stems from both physical and numerical differences between the codes, which we will now alleviate by adjusting the CONCEPT simulation.

Cell-centred vs. cell-vertex discretisation

The initial particle distribution, stored within the generated initial snapshot file ic/tutorial-9_a=0.02, should ideally be close to homogeneous and isotropic. We can check this qualitatively by plotting the distribution in 3D. For this, let’s use the render3D utility:

./concept \
    -u render3D ic/tutorial-9_a=0.02 \
    -p param/tutorial-9 \
    -c "render3D_options = { \
        'interpolation': 0, \
        'color': 'black', \
        'background': 'white', \
        'resolution': 2000}"

This produces an image file in the ic directory. Zooming in on one of the corners of the box, it should be clear that while we do have a close to homogeneous system, it is far from isotropic.

Note

Adding the render3D_options parameter above tunes the render so that it more clearly shows what we’re after. Indeed, a more interesting render is obtained by not adding this parameter.

The strong anisotropy of the particles is inherited from the grid on which they are (pre-)initialised. This grid-like structure of the particle distribution makes the system sensitive towards other grids used during the simulation, e.g. the potential grid. In CONCEPT, a default choice of using ‘cell-centred’ grid discretisation is made, whereas GADGET-2 uses ‘cell-vertex’ grid discretisation. In effect this means that the potential grids of the two codes are shifted by half a grid cell relative to each other, in all three dimensions (see the cell_centered parameter for more information). This detail oughtn’t matter much, but the grid structure imprinted on the particle distribution, together with the fact that our simulations are rather small (\(N = 64^3\)), makes for an exaggerated effect.

We can force CONCEPT to use cell-vertex discretisation by specifying

cell_centered = False

Introduce this new parameter — either by editing the parameter file or supplying it with -c — and rerun the CONCEPT simulation. After also updating the comparison plot, you should find that the relative error is now much smaller, though still a few percent.

Note

Though we now have better agreement, we should not be tempted to conclude that cell-vertex discretisation is superior to its cell-centred cousin. The improvement came from using the same discretisation.

The observed difference in results for the two discretisation schemes would be smaller had the initial conditions not contained the strong anisotropies. CONCEPT does allow for the creation of more isotropic initial condition by initially placing the particles on either a body-centered cubic (bcc) lattice or a face-centered cubic (fcc) lattice, rather than the standard simple cubic (sc) lattice. See the initial_conditions parameter for more information.

Though undesirable when comparing against GADGET-2, CONCEPT additionally effectively allows for simultaneous usage of both discretisation schemes by carrying out both of them followed by an interlacing step to combine the results. This technique effectively lowers the discretisation scale, reducing the numerical artefacts. For more information about potential interlacing, see the potential_options parameter.

Radiation in the background

A key difference between CONCEPT and GADGET-2 is that the former — unlike the latter — tries to incorporate all species of a given cosmology, as demonstrated in the previous section. As such, CONCEPT uses a full background evolution obtained from the CLASS code. By default, this CLASS background includes matter, \(\Lambda\) and radiation, whereas the background of GADGET-2 only consists of matter and \(\Lambda\).

One cannot remove radiation from the CLASS background, as a cosmology without radiation is not supported by CLASS. Instead, we can deactivate the CLASS background entirely, which will turn on an internal background solver of CONCEPT, which like that of GADGET-2 only includes matter and \(\Lambda\). To do so, add

enable_class_background = False

to the parameters (while also still keeping cell_centered = False) and rerun the CONCEPT simulation.

Updating the comparison plot, the codes should now agree to within a percent.

Increased short-range precision

The gravitational force of the P³M method in CONCEPT and the TreePM method of GADGET-2 is split into a long-range and a short-range part at a scale \(x_{\text{s}}\), which in both CONCEPT and GADGET-2 has a default value of \(x_{\text{s}} = 1.25 \Delta x\), with \(\Delta x\) the cell size of the potential grid. The short-range force between pairs of particles separated by a distance somewhat larger than \(x_{\text{s}}\) falls off exponentially with the distance. Particle pairs with a separation \(|\boldsymbol{x}_i - \boldsymbol{x}_j| > x_{\text{r}}\) may then be neglected, provided \(x_{\text{r}} \gg x_{\text{s}}\).

In both CONCEPT and GADGET-2, the range of the short-range force has a default value of \(x_{\text{r}} = 4.5x_{\text{s}}\). Due to the grouping of particles by the tree in GADGET-2, the short-range forces between some particles separated by distances somewhat larger than \(x_{\text{r}}\) are being taken into account as well. We can then hope to obtain still better agreement with GADGET-2 by increasing the value of \(x_{\text{r}}\) in CONCEPT slightly. Try rerunning CONCEPT with \(x_{\text{r}} = 5.5x_{\text{s}}\) with the added parameter

shortrange_params = {'range': '5.5*scale'}

and updating the comparison plot. You should find that the two codes now agree well within 1%.

For still better agreement, you may try similarly upgrading \(x_{\text{r}}\) to \(5.5x_{\text{s}}\) in GADGET-2 (either by including the updated shortrange_params when running the CONCEPT gadget utility or by manually setting RCUT in output/tutorial-9/Gadget2/Makefile) and then re-compiling and -running GADGET-2 and computing the power spectra from the new snapshots. Improved agreement between the two codes is also achieved by increasing the simulation size.

Final notes

CONCEPT strives to be consistent with general relativistic perturbation theory. It does so by using the full background from CLASS, and adding in gravitational effects from linear perturbations of other species (when enabled), ensuring that the simulation stays in the so-called N-body gauge. All this fanciness is hardly needed when running without massive neutrinos or other exotic species, which GADGET-2 does not support.

GADGET-2 only allows for the simulation of matter (though besides the standard cold dark matter particles, it further supports (smoothed-particle) hydrodynamical baryons), and always uses a background containing just matter and \(\Lambda\). Sticking to even the simplest ΛCDM cosmologies (which do include photons), neglecting the gravitational tug from radiation is unacceptable for precision simulations. To resolve this problem, the usual trick used with Newtonian N-body codes is to generate initial conditions at \(a = a_{\text{begin}}\) using the results of general relativistic perturbations theory (which also accounts for radiation) at \(a = 1\), but scaled back in time using the Newtonian growth factor. Thus, only the \(a = 1\) results of such back-scaled simulations are really physical. The realization_options parameter in the parameter file used for this section specifies that back-scaling should be used when generating the initial conditions, and also that this should be done within the synchronous gauge (as opposed to the default N-body gauge), as is typically used with back-scaling.

In the above, we eventually disabled the CLASS background within CONCEPT simulations by setting enable_class_background = False, which then instead makes use of a simplified matter + Λ background, similar to what is used within GADGET-2. Note however that we did this after creating the initial conditions, meaning that the back-scaling was done in a background that includes radiation, causing a slight inconsistency. For the comparison between CONCEPT and GADGET-2, this is not relevant, but for physically accurate results, the same value of enable_class_background should be used when generating the initial conditions as used within the simulation itself.

When using back-scaling, the particle velocities are obtained from the same (\(\delta\)) transfer function as the positions, using a Newtonian approximation. The velocities are thus not the actual synchronous gauge velocities, but are in fact closer to the velocities within the N-body gauge, as obtained from the \(\theta\) transfer function and no back-scaling. As the synchronous and N-body gauge \(\delta\) transfer functions are quite similar, the back-scaled synchronous gauge initial conditions — used within this section — are then in fact very similar to the non-back-scaled N-body gauge initial conditions used within CONCEPT by default. The specifications within realization_options in the parameter file are thus not very crucial for the results, but mostly put there to comply with the spirit of Newtonian N-body simulations and GADGET-2. We note again that the real difference between the N-body gauge approach and the back-scaled, synchronous gauge approach is to be found outside basic ΛCDM cosmologies and matter-only simulations, where the Newtonian approximations used for the back-scaling break down.