Running and post-processing¶

GPU v2.0 users

This page is oriented toward the CPU/v1.x documentation. For the GPU release, use the gevolution 2.0 running guide and output guide.

This section walks through a minimal end-to-end run and then shows how to inspect outputs. The examples assume you run from the project root and that your output directory is defined in your settings (typically in settings.ini). This example shows you the most basic way to interact with the output.

Running a basic simulation¶

Here is how to run a simple simulation:

Configure the makefile and build the code: make clean; make. See the documentation.
Configure your run in settings.ini (box size, resolution, output cadence, lightcone options). See the settings file guide and the documentation
Run the code with MPI and confirm the output directory is created. The output path is set in the settings file. The code can be run by, for example,
```
    mpirun -np 4 ./gevolution -n 2 -m 2 -s settings.ini 
```
Post-process with the tools below.

Tip: Either differentiate different simulations by varying their paths, or the generic file base parameter in settings file.

Processing the output¶

Below are inlined Jupyter notebook-style examples. You can paste each cell into a notebook to explore the outputs interactively.

Notebook setup cell

We start by loading some useful modules and setting the path.

import matplotlib.pyplot as plt, numpy as np

# path to simulation
path = './outputs/'

For all of the below, we will assume that generic file base, snapshot file base and Pk file base are all left to their default values, see the settings file guide.

Power spectra¶

Power spectra are calculated on the fly and typically saved \(\sim 50\) times for one simulation, due to their small file sizes. The power spectra are dimensionless. Thus the units have to be converted for comparison with e.g. CLASS. The example below looks for a spectrum file and plots \(P_m(k)\).

file = 'pk000_delta.dat'
pk = np.loadtxt(path+file)
! head -n 3 {path+file} # see the redshift + column names
plt.plot(pk[:,0],pk[:,1]/pk[:,0]**3) # CLASS format
plt.xlabel(r'$h$ Mpc$^-1$'); plt.ylabel(r'$(h^{-1}$Mpc)^3$')

Tip: Getting the string for the snapshot/pk/etc number can be done with the function
    getnum = lambda num: '0'*int(2-int(np.log10(max(num,1))))+str(num)

CDM snapshot¶

Particle snapshots are stored in Gadget2 binary format. For large runs, they should be requested as snapshot outputs = multi-Gadget2, that makes each MPI processor outputs its domain into a file snapshot<number>_<species>.<proc_number>. The cell below shows a common Gadget2-reading workflow for inspecting positions and velocities. We make use of the open-source module pygadgetreader.

a = 1 # assumed redshift z=0
import pygadgetreader as gr
file = 'snapshot000_cdm.0'
pos = gr.readsnap(path+file,'pos',1)/1e3 # now in Mpc/h
vel = gr.readsnap(path+file,'vel',1)/a**0.5 # in physical km/s
ids = gr.readsnap(path+file,'id',1) 
# plot the projected particle distribution in the XY plane (at most 10k points for speed)
mx = int(1e4); stride = max(float(len(pos))/float(mx),1)
plt.scatter(pos[::stride,0],pos[::stride,1])
plt.xlabel(r'$Mpc/h$'); plt.ylabel(r'$Mpc/h$')

Field snapshots¶

Field outputs (e.g. density, potential, metric perturbations) are stored on grids. The following cell demonstrates loading a 3D field and plotting a 2D slice. Due to difference between row-major and column-major ordering, we need to transpose the fields for comparison with the particle snapshots above.

import h5py

def geth5arr(filename, ind=None):
    with h5py.File(filename, 'r') as infile:
        keys = list(infile.keys())
        if not keys:
            raise ValueError("HDF5 file has no datasets.")
        key = keys[0]
        data_raw = np.array(infile[key])
        if ind is None:
            data = data_raw.transpose(2, 0, 1)
        else:
            data = data_raw[:, :, :, ind].transpose(2, 1, 0)
    return data

filename = 'snapshot000_T00.h5'
boxsize = 1000
T00 = geth5arr(path+file)[:,:,0] # scalar, so ind=None
plt.imshow(T00/np.mean(T00)-1.,extent=[0,boxsize,0,boxsize]) # matter overdensity plot
plt.colorbar(label=r'$\delta_m$')
plt.xlabel(r'$h^{-1}$Mpc'); plt.ylabel(r'$h^{-1}$Mpc')

Lightcone data¶

Lightcone outputs collect data on the past lightcone. They are stored in the binary .map format, divided in shells with some overlap for ease of interpolation. They can be integrated and converted to fits maps in healpix format using the lcmap software provided as part of gevolution. lcmap outputs the gravitational lensing potential kappa, the Integrated-Sachs-Wolfe (ISW) contribution isw, the Sachs-Wolfe (SW) contribution phi and the Shapiro time delay shapiro.

First make and run lcmap. We assume a pencil beam lightcone number 1 until a distance 1500 \(h^{-1}\)Mpc, and a fullsky lightcone number 0 until 200 \(h^{-1}\)Mpc.

make clean; make lcmap
./lcmap -d 1500 -l 0,1 -s ./outputs/_settings_used.ini

where there should not be a space between 0, and 1.

Tip: A pencil-beam lightcone can extend past the size of the box volume by use of the periodic boundary without self-intersections by choosing the angle and opening carefully. For narrow beams, a supplementary small full-sky lightcone close to the origin is recommended, to avoid issues with sampling the fields close to origin. lcmap will automatically switch between the lightcones in the integration.

We can now read the fits file as a healpix map, using healpy and astropy. Let's plot the ISW map:

import healpy as hp
from astropy.io import fits
file = 'lightcone_isw.fits'
pix = np.array(fits.open(path+file)[1].data)
hp.gnomview(p,rot=(90,90,0), xsize=800,title='',cbar=True) # adjust parameters depending on lightcone