class utility

The CONCEPT ‘class’ utility enables easy access to linear perturbations (and background quantities) outside of CONCEPT, e.g. for use with other simulation codes. In a nutshell, a cosmology is specified to CONCEPT using the standard parameter specification, with which CONCEPT calls CLASS in order to obtain the required perturbation data. All perturbations are then processed; interpolating them onto a common \((a, k)\) grid, and typically also transforming them to N-body gauge. The resulting data is saved to an HDF5 file, ready to be used for e.g. initial condition generation or gravitational corrections on matter particles from linear species such as neutrinos. The pkdgrav3 code is able to fully utilise such HDF5 files generated by the CONCEPT class utility, enabling simulations fully consistent with general relativistic perturbation theory.

Tip

If you are unfamiliar with the basic framework of the N-body gauge for cosmological simulations, you might want to check out the paper on ‘Fully relativistic treatment of light neutrinos in 𝘕-body simulations’. It is also advisable to have gone through the CONCEPT tutorial at least up to (and including) the first subsection of beyond matter-only simulations.

Entries on this page:

Basic usage 

A basic invocation of the CONCEPT class utility looks like

./concept \
    -u class \
        --kmin <kmin> --kmax <kmax> \
        <perturbations>

where <kmin> and <kmax> specify the \(k\) range while <perturbations> is a string of comma-separated perturbation specifications. As an example, the below processes \(\delta\) and \(\theta\) perturbations for ‘total matter’ (combined baryons and cold dark matter), as well as \(\delta\) perturbations for photons and massless neutrinos, in the interval \(10^{-3}\,\text{Mpc}^{-1} \leq k \leq 1 \text{Mpc}^{-1}\):

./concept \
    -u class \
        --kmin 1e-3/Mpc --kmax 1/Mpc \
        "b+cdm:1, g:0, ur:0"

Each specification within <perturbations> consists of a CLASS species and a “Boltzmann order” separated by a colon. Commonly used CLASS species include

b: Baryons.
cdm: Cold dark matter.
g: Photons.
ur: Massless neutrinos.
ncdm[0], ncdm[1], ncdm[2], … : Massive neutrinos.
metric: Fictitious species providing relativistic corrections (not a standard CLASS species, provided by CONCEPT).

As in the above example, two (or more) species may be joined using +, leading to combined, weighted perturbations.

The available Boltzmann orders are

0: Process \(\delta\) (energy density contrast) perturbations.
1: Process \(\delta\) and \(\theta\) (velocity divergence) perturbations.
2: Process \(\delta\), \(\theta\), \(\delta P\) (pressure) and \(\sigma\) (shear stress) perturbations.

The Boltzmann order defaults to 0 if left out.

Running the CONCEPT class utility as above will produce a file named class_processed.hdf5 within the output directory, containing the specified perturbations as well as various background quantities. The contents of these HDF5 files can be explored using tools such as ViTables, but is also described below.

For selection of the times (scale factor \(a\) values), see the --times command-line option.

For selection of the \(k\) modes, see the --modes command-line option.

Command-line options 

Below the possible command-line options to the CONCEPT class utility are described, i.e. additional options to the concept script allowed once -u class is supplied. For the <perturbations> option, see basic usage above. See this page for the command-line options to the concept script in general, and specifically this section for how standard command-line options are used together with utility specific ones.

Help: `-h`, `--help`

Displays a short description of each command-line option to the class utility and exits:

./concept -u class -h

Minimum Fourier mode: `--kmin`

Sets the minimum needed Fourier \(k\) mode for the perturbations. For e.g. \(k \geq 5\times 10^{-4}\, \text{Mpc}^{-1}\):

./concept -u class --kmin 5e-4/Mpc

Note

A small tolerance is automatically applied so that the actual lowest \(k\) mode will be slightly below what is specified.

Maximum Fourier mode: `--kmax`

Sets the maximum needed Fourier \(k\) mode for the perturbations. For e.g. \(k \leq 2\, \text{Mpc}^{-1}\):

./concept -u class --kmax 2/Mpc

Note

A small tolerance is automatically applied so that the actual largest \(k\) mode will be slightly above what is specified.

Tip

Try not to specify a maximum \(k\) much larger than what you need, as the CLASS computation time for each mode grows rapidly with \(k\).

Fourier modes: `--modes`

When specified as an integer, this sets the total number of Fourier \(k\) modes for the perturbations. E.g.

./concept -u class --modes 128

The placement of the \(k\) modes at which to tabulate the perturbations are chosen based on the minimum and maximum mode together with the class_modes_per_decade parameter. When the --modes option is not specified, class_modes_per_decade is used as is, leading to some number of total modes. When the --modes option is specified, the numbers of modes per decade within class_modes_per_decade are all uniformly scaled in order to yield the requested number of modes.

Tip

If you simply want some definite number (say \(128\)) of logarithmically equidistant modes, you can do e.g.

./concept -u class --modes 128 -c "class_modes_per_decade = 1"

Tip

If the processed perturbations are intended for use with pkdgrav3, note that a maximum number of \(256\) modes is allowed.

Other allowed specifications for --modes are:

Explicit \(k\) values at which to perform the tabulation, specified as any valid Python expression, e.g.
```
./concept -u class --modes "array([1e-2, 1e-1, 1])/Mpc"
```
Path to text file containing values of \(k\). It is assumed that the data within this file is expressed in the same units as used by the current CONCEPT run.
Path to HDF5 file produced using the CONCEPT class utility. The \(k\) values to use will be taken from the tabulated perturbations within this file.

In all cases, the minimum and maximum modes are respected, with \(k\) values outside this range being dropped.

Times: `--times`

When specified as an integer, this sets the total number of times (scale factor \(a\) values) at which to tabulate the perturbations. E.g.

./concept -u class --times 512

When not specified, this defaults to 1024 times. The placement of the tabulation times follows that which is used by CLASS, providing denser sampling around epochs of greater change.

Note

If the raw perturbations from the CLASS computation are tabulated at fewer points in time than requested, no additional times will be added. It is thus possible that the processed perturbations are tabulated at fewer times than requested through --times (though this should not happen for reasonable values).

The \(a\) values will be within the interval \(a_{\text{begin}} \leq a \leq 1\), with the a_begin parameter specifying the lower boundary.

Note

A small tolerance is automatically applied so that the actual smallest \(a\) value will be slightly below what is specified through a_begin.

Tip

If the processed perturbations are intended for use with pkdgrav3, note that a maximum number of \(1024\) times is allowed.

Other allowed specifications for --times are:

Explicit \(a\) values at which to perform the tabulation, specified as any valid Python expression, e.g.
```
./concept -u class --times "[1/(1 + z) for z in linspace(0, 99, 512)]"
```
Path to text file containing values of \(a\).
Path to HDF5 file produced using the CONCEPT class utility. The \(a\) values to use will be taken from the tabulated perturbations within this file.

Gauge: `--gauge`

Sets the gauge of the processed perturbations. E.g.

./concept -u class --gauge N-body

The available gauges are:

N-body: The N-body gauge (default).
synchronous: The synchronous gauge.
Newtonian: The conformal Newtonian (also know as the longitudinal) gauge.

When the N-body gauge is selected, CLASS is really run in synchronous gauge, with the perturbations transformed to N-body gauge by CONCEPT.

File contents 

The result of a CONCEPT class utility computation is a file called class_processed.hdf5, containing the requested perturbations as well as background quantities. As mentioned earlier, the contents of HDF5 files are conveniently explored using a graphical tool like ViTables. Details of the contents of class_processed.hdf5 are described below.

Data layout

The data layout within class_processed.hdf5 is as follows:

background: Group containing data sets for various background quantities, tabulated at a common 1D time grid:
- a: Scale factor \(a\) values at which all background quantities are tabulated.
- z: Redshift \(z\) values at which all background quantities are tabulated.
- t: Cosmic time \(t\) values at which all background quantities are tabulated.
- tau: Conformal time \(\tau\). values at which all background quantities are tabulated.
- H: The Hubble parameter \(H\).
- rho_<species> with <species> some CLASS species: The background energy density \(\bar{\rho}\) of the given species. The critical ‘crit’ and total ‘tot’ energy densities are present as well.
- p_<species> with <species> some CLASS species: The background pressure \(\bar{P}\) of the given species.
- D: First-order (linear) growth factor \(D = D^{(1)}\), with the convention \(D^{(1)}(a=1) = 1\).
- f: First-order (linear) growth rate \(f = f^{(1)} \equiv \mathrm{d}\ln D^{(1)} / \mathrm{d}\ln a\).
- D2: Second-order growth factor \(D^{(2)}\), with the convention \(D^{(2)} > 0\).
- f2: Second-order growth rate \(f^{(2)} \equiv \mathrm{d}\ln D^{(2)} / \mathrm{d}\ln a\).
- D3a: Third-order growth factor \(D^{(3\text{a})}\), with the convention \(D^{(3\text{a})} > 0\).
- f3a: Third-order growth rate \(f^{(3\text{a})} \equiv \mathrm{d}\ln D^{(3\text{a})} / \mathrm{d}\ln a\).
- D3b: Third-order growth factor \(D^{(3\text{b})}\), with the convention \(D^{(3\text{b})} > 0\).
- f3a: Third-order growth rate \(f^{(3\text{b})} \equiv \mathrm{d}\ln D^{(3\text{b})} / \mathrm{d}\ln a\).
- D3c: Third-order growth factor \(D^{(3\text{c})}\), with the convention \(D^{(3\text{c})} > 0\).
- f3c: Third-order growth rate \(f^{(3\text{c})} \equiv \mathrm{d}\ln D^{(3\text{c})} / \mathrm{d}\ln a\).
In addition to the above datasets, this group also stores some attributes, namely the reduced Hubble constant h (given by \(h \equiv H_0/(100\, \text{km}\, \text{s}^{-1}\, \text{Mpc}^{-1})\) as well as density parameters Omega_<species> for the various species.
perturbations: Group containing the specified perturbations, tabulated at a common 2D \((a, k)\) grid:
- a: Data set of scale factor \(a\) values at which all perturbations are tabulated.
- k: Data set of Fourier modes \(k\) at which all perturbations are tabulated.
- delta_<species>: The \(\delta\) perturbations of the given species.
- theta_<species>: The \(\theta\) perturbations of the given species.
- deltaP_<species>: The \(\delta P\) perturbations of the given species.
- sigma_<species>: The \(\sigma\) perturbations of the given species.
In addition to the above datasets, this group also stores a gauge attribute, specifying the gauge of the perturbations.
class_params: Group containing no data sets, but with all CLASS parameters used for the CLASS computation stored as attributes. That is, the attributes within this group contain all information needed to do a CLASS run identical to the one carried out by the CONCEPT class utility.

Note

One of these CLASS parameters is gauge. When running the CONCEPT class utility in N-body gauge (the default), CLASS is really run in synchronous gauge (with the results transformed to N-body gauge by CONCEPT), and so the gauge attribute inside the class_params group will be set to 'synchronous'. For the actual gauge in which the perturbations are expressed, see the gauge attribute on the perturbations group.
units: Group containing no data sets, but with attributes storing the units used for the data throughout the file:
- unit length: The length unit used. By default this is Mpc.
- unit time: The time unit used. By default this is Gyr.
- unit mass: The time unit used. By default this is 10**10*m_sun (\(10^{10}\, m_{\odot}\)).

Extra output

The above lists of background and perturbation quantities are what end up in the processed HDF5 file by default, though the raw CLASS computation contains other quantities as well. You may request these other quantities by listing them in the class_extra_background parameter and the class_extra_perturbations parameter.

As an example, the below requests the comoving distance as an extra background quantity, as well as the conformal Newtonian metric potentials \(\phi\) and \(\psi\) as extra perturbations:

./concept \
    -u class \
        --kmin <kmin> --kmax <kmax> \
        <perturbations> \    -c "class_extra_background = 'comov. dist.'" \    -c "class_extra_perturbations = {'ϕ', 'ψ'}"

Note

While the various perturbations in general depend on the chosen gauge, this is not so for the gauge specific metric perturbations ϕ, ψ, hʹ and H_Tʹ, which are always given in their definite gauges. Also, hʹ is not available when running in conformal Newtonian gauge.

Caution

The units of the various quantities within class_processed.hdf5 are described below. If adding an extra CLASS quantity via class_extra_background or class_extra_perturbations that is unknown to CONCEPT, unit convertion will not be carried out, and the quantity will be provided exactly as CONCEPT got it from CLASS. If this happens, a warning will be emitted.

Units

As written above, the units group contains the base units used for all data within the file, with the default base units being \(\text{Mpc}\), \(\text{Gpc}\) and \(10^{10}\, m_{\odot}\). Given these default base units, the table below lists the full units in which the various data sets are expressed:

Quantity	Unit (default)
Scale factor `a` and redshift `z`	\([1]\)
Cosmic time `t` and conformal time `tau`	\([\text{Gyr}]\)
Fourier mode `k`	\([\text{Mpc}^{-1}]\)
Hubble parameter `H`	\([\text{Gyr}^{-1}]\)
Growth factors `D`, `D2`, `D3a`, `D3b`, `D3c` and rates `f`, `f2`, `f3a`, `f3b`, `f3c`	\([1]\)
Background energy density `rho`	\([10^{10}\, m_{\odot}\, \text{Mpc}^{-3}]\)
Energy density contrast `delta`	\([1]\)
Velocity divergence `theta`	\([\text{Gyr}^{-1}]\)
Pressure background `p` and perturbation `deltaP`	\([10^{10}\, m_{\odot}\, \text{Mpc}^{-1}\,\text{Gyr}^{-2}]\)
Shear stress `sigma`	\([\text{Mpc}^2\,\text{Gyr}^{-2}]\)
Metric perturbations `phi` and `psi`	\([\text{Mpc}^2\,\text{Gyr}^{-2}]\)
Time derivatives of metric perturbations `h_prime` and `H_T_prime`	\([\text{Gyr}^{-1}]\)

Tip

If you need the Hubble parameter H in its canonial units of \(\text{km}\, \text{s}^{-1}\, \text{Mpc}^{-1}\), you can convert it by multiplying the values as follows:

\[\begin{split}\require{upgreek} H \rightarrow 1\times H &\equiv \frac{1495978707}{487000\uppi}\, \frac{\text{km}\, \text{s}^{-1}\, \text{Mpc}^{-1}}{\text{Gyr}^{-1}} H \\ &= 977.7922216807892\, \frac{\text{km}\, \text{s}^{-1}\, \text{Mpc}^{-1}}{\text{Gyr}^{-1}} H\end{split}\]

Tip

The background energy densities rho are really expressed in units of mass densities. If you need these in energy density units, multiply by the speed of light squared:

\[\begin{split}\require{upgreek} \rho \rightarrow \rho c^2\,, \quad c &\equiv \frac{1999985302\uppi}{20492859}\, \text{Mpc}\, \text{Gyr}^{-1} \\ &= 306.60139378555056\, \text{Mpc}\, \text{Gyr}^{-1}\end{split}\]

Note

The base unit system in which to express the data within class_processed.hdf5 may be changed through the unit_length, unit_time and unit_mass parameters.

Complete example (massive neutrinos)

Here follows a complete code example for a CONCEPT class utility computation of a cosmology involving massive neutrinos.

The cosmology for which we want to compute the linear perturbations is specified by the following table:

Parameter	Value
\(H_0\)	\(67\, \text{km}\, \text{s}^{-1}\, \text{Mpc}^{-1}\)
\(\Omega_{\text{b}}\)	\(0.049\)
\(\Omega_{\text{cdm}} + \Omega_{\nu}\)	\(0.27\)
\(T_{\gamma}\)	\(2.7255\, \text{K}\)
\(N_{\text{eff}}\)	\(3.046\)
\(m_{\nu_{1}}\)	\(0\, \text{eV}\)
\(m_{\nu_{2}}\)	\(8.7\times 10^{-3}\, \text{eV}\)
\(m_{\nu_{3}}\)	\(5.0\times 10^{-2}\, \text{eV}\)

The total neutrino density parameter \(\Omega_{\nu} = \Omega_{\nu_1} + \Omega_{\nu_2} + \Omega_{\nu_3}\) is implicitly set through the effective number of neutrino species \(N_{\text{eff}}\) together with the neutrino masses \(m_1\), \(m_2\), \(m_3\). The neutrino temperature \(T_{\nu}\) (equal for all three species) is related to the photon temperature \(T_{\gamma}\) by

\[T_{\nu} = \biggl(\frac{4}{11}\biggr)^{\frac{1}{3}} \biggl(\frac{N_{\text{eff}}}{3}\biggr)^{\frac{1}{4}} T_{\gamma} \,.\]

Note

We do not need to specify any parameters having to do with the primordial spectrum (e.g. \(A_{\text{s}}\) and \(n_{\text{s}}\)) as these do not enter the linear CLASS computation.

The below parameter file implements the above neutrino cosmology, though written in a way that allows one to easily change the number of neutrino species:

class-param

# Input/output
output_dirs = param.dir

# Cosmology
H0 = 67*km/(s*Mpc)
Ωb = 0.049
Ωcdm = 0.27 - Ων
a_begin = 1/(1 + 99)
_mν = [0, 8.7e-3, 5.0e-2]  # neutrino masses in eV
_N_eff = 3.046
class_params = {
    # Photon temperature
    'T_cmb': 2.7255,
    # Neutrino hierarchy
    'N_ur'    : 0,
    'N_ncdm'  : len(set(_mν)),
    'deg_ncdm': [_mν.count(mν) for mν in sorted(set(_mν))],
    'm_ncdm'  : [mν if mν > 0 else 1e-100 for mν in sorted(set(_mν))],
    'T_ncdm'  : [(4/11)**(1/3)*(_N_eff/len(_mν))**(1/4)]*len(set(_mν)),
    # Photon precision parameters
    'l_max_g'                          : 1000,
    'l_max_pol_g'                      : 1000,
    'radiation_streaming_approximation': 3,
    # Massive neutrino precision parameters
    'l_max_ncdm'              : 200,
    'Number of momentum bins' : [50]*len(set(_mν)),
    'Quadrature strategy'     : [2]*len(set(_mν)),
    'ncdm_fluid_approximation': 3,
    # General precision parameters
    'evolver'                     : 0,
    'recfast_Nz0'                 : 1e+5,
    'tol_thermo_integration'      : 1e-6,
    'perturb_integration_stepsize': 0.25,
    'perturb_sampling_stepsize'   : 0.01,
}

The key points of this parameter file are:

Setting output_dirs to param.dir ensures that the resulting HDF5 file with the perturbations will be placed in the same directory as the parameter file.
The amount of cold dark matter \(\Omega_{\text{cdm}}\) is set through Ωcdm = 0.27 - Ων, though Ων (\(\Omega_{\nu}\)) is not explicitly set. As stated previously, \(\Omega_{\nu}\) is implicitly defined through the other neutrino parameters, from which CONCEPT automatically sets Ων.
By specifying a_begin = 1/(1 + 99) we declare that we are not interested in perturbations prior to \(z = 99\).
The neutrino masses _mν and effective number of neutrino species _N_eff are both stored in variables starting with an underscore ‘_’. This is the canonical naming convention within CONCEPT parameter files for helper variables that are not supposed to be CONCEPT parameters.
The class_params parameter holds parameters fed to CLASS.
- Note that 'H0', 'Omega_b' and 'Omega_cdm' is automatically added from H0, Ωb and Ωcdm.
- The parameters are set up to take advantage of the case where _mν holds several identical neutrino masses (as in e.g. _mν = [0, 3e-2, 3e-2]). In this case, only neutrinos with unique masses (set(_mν)) will be computed (saving computation time), with their degeneracy recorded in the 'deg_ncdm' CLASS parameter.
- We use the ncdm species for all neutrinos, including the massless one. We need to explicitly set the number of massless neutrinos, N_ur, to zero, as this is non-zero by default. To avoid possible issues within the CLASS computation from having a mass identically equal to zero using ncdm, we replace massless entries in 'm_ncdm' with \(10^{-100}\).
- In standard CLASS, parameters with multiple values (e.g. 'm_ncdm') needs to be specified as a comma-separated string. In CONCEPT parameter files, using lists ‘[...]’ is possible as well.
- The many precision parameters improve the correctness of various perturbation results significantly over the default, but also makes the CLASS computation much more expensive.

Tip

If you would like to better understand a given CLASS parameter, you can look for it in the explanatory.ini example CLASS parameter file, which documents many of the possible CLASS parameters.

With the above parameter file saved to a file called e.g. class-param, we can run the CONCEPT class utility on it using something like

./concept \
    -u class \
        --kmin "1e-3*0.67/Mpc" --kmax "1*0.67/Mpc" --modes 128 \
        --times 512 \
        "b+cdm:1, g, ncdm[0], ncdm[1], ncdm[2]:2, metric" \
    -p /path/to/class-param \
    -n 4

Tip

If working on a cluster, you can automatically submit the computation as a job, just as when running full simulations.

The above command requests \(128\) modes between \(k = 10^{-3}\, h\, \text{Mpc}^{-1}\) and \(k = 1\, h\, \text{Mpc}^{-1}\), with \(h = H_0/(100\, \text{km}\, \text{s}^{-1}\, \text{Mpc}^{-1}) = 0.67\), all tabulated at the same \(512\) times. Both \(\delta\) and \(\theta\) perturbations of total matter are requested, allowing for generation of N-body initial conditions. Furthermore, \(\delta\) perturbations in all linear species — photons, the three neutrino species and the metric correction species — are requested, which might be used to supply linear gravitational kicks to matter particles within an N-body simulation. Lastly, \(\theta\), \(\delta P\) and \(\sigma\) perturbations for the heaviest neutrino species are also requested. The perturbations will all be in N-body gauge. Using 4 processes as in the above example command, the entire computation will take around 6 minutes on modern hardware.

With the class_processed.hdf5 file generated from the above command and parameter file, the below Python script demonstrates how one can extract and use (some of) the data within this file:

plot.py

import re
import numpy as np
import matplotlib; matplotlib.use('agg')
import matplotlib.pyplot as plt
import scipy.interpolate
import h5py

# Read in data from the HDF5 file
with h5py.File('class_processed.hdf5', mode='r') as f:
    # Read in all background energy densities
    background = f['background']
    print('Available background quantities:', list(background.keys()))
    a_bg = background['a'][:]
    ρ_bg = {}
    for key, val in background.items():
        match = re.search(r'^rho_(.+)$', key)
        if match:
            ρ_bg[match.group(1)] = val[:]
    h = background.attrs['h']
    # Read in all δ perturbations
    perturbations = f['perturbations']
    print('Available perturbations:', list(perturbations.keys()))
    a_pt = perturbations['a'][:]
    k_pt = perturbations['k'][:]
    δ_pt = {}
    for key, val in perturbations.items():
        match = re.search(r'^delta_(.+)$', key)
        if match:
            δ_pt[match.group(1)] = val[...]
    gauge = perturbations.attrs['gauge']
    # Get unit system
    units = dict(f['units'].attrs)
    print('Units:', units)

# Construct matter power spectrum interpolator
A_s = 2.1e-9
n_s = 0.96
α_s = 0
k_pivot = 0.05  # [Mpc⁻¹]
ζ = lambda k: (
    np.pi*np.sqrt(2*A_s)
    *k**(-3/2)*(k/k_pivot)**((n_s - 1)/2)
    *np.exp(α_s/4*np.log(k/k_pivot)**2)
)
interpolator = scipy.interpolate.RegularGridInterpolator(
    (np.log(k_pt), np.log(a_pt)),
    np.log((ζ(k_pt)*δ_pt['b+cdm'])**2).T,
    method='cubic',
)
P_m = lambda z, logk=np.log(k_pt), interpolator=interpolator: (
    np.exp(interpolator((logk, np.log(1/(1 + z)))))
)

# Construct δρ interpolators for all species
δρ = {}
for species, δ_species in δ_pt.items():
    ρ_species = np.exp(scipy.interpolate.interp1d(
        np.log(a_bg),
        np.log(ρ_bg[species]),
        kind='cubic',
    )(np.log(a_pt)))
    interpolator = scipy.interpolate.RegularGridInterpolator(
        (np.log(k_pt), np.log(a_pt)),
        (δ_species*ρ_species.reshape(-1, 1)).T,
        method='linear',
    )
    δρ[species] = lambda z, logk=np.log(k_pt), interpolator=interpolator: (
        interpolator((logk, np.log(1/(1 + z))))
    )

# Plot energy density history
fig, axes = plt.subplots(2, 2, figsize=(9, 6.75))
ax = axes[0, 0]
a_min = 1e-3
mask = (a_bg >= a_min)
for species, ρ in ρ_bg.items():
    match = re.search(r'^(b|cdm|g|ncdm\[\d+\]|lambda)$', species)
    if match:
        ax.loglog(a_bg[mask], ρ[mask]/ρ_bg['crit'][mask], label=species)
ax.set_xlim(a_min, 1)
ax.set_ylim(bottom=1e-5)
ax.legend(fontsize=8)
ax.set_xlabel(r'$a$')
ax.set_ylabel(r'$\bar{{\rho}}_{\alpha}/\rho_{\mathrm{crit}}$')
ax.set_title(f'Energy density history')

# Plot matter power spectrum
ax = axes[0, 1]
z_values = [0, 0.5, 1, 2, 5, 10]
for z in z_values:
    ax.loglog(k_pt/h, h**3*P_m(z), label=rf'$z = {z}$')
ax.set_xlim(k_pt[0]/h, k_pt[-1]/h)
ax.legend(fontsize=8)
unit_length = units['unit length']
ax.set_xlabel(rf'$k\, \left[h\, \mathrm{{{unit_length}}}^{{-1}}\right]$')
ax.set_ylabel(
    rf'$P_{{\mathrm{{b}}+\mathrm{{cdm}}}}\, '
    rf'\left[(h^{{-1}}\, \mathrm{{{unit_length}}})^3\right]$'
)
ax.yaxis.tick_right()
ax.yaxis.set_label_position('right')
ax.set_title(f'Matter power spectrum ({gauge} gauge)')

# Plot δρ perturbations of all species
z_values = [50, 1]
unit_mass = (
    units['unit mass']
    .replace(r'**(10)', r'^{10}')
    .replace(r'*m_sun', r'\, m_{\odot}')
)
for ax, z in zip(axes[1, :], z_values):
    for i, (species, δρ_species) in enumerate(δρ.items()):
        for sign, linestyle, label in zip([-1, +1], ['-', ':'], [species, None]):
            ax.loglog(k_pt/h, sign*δρ_species(z), f'C{i}{linestyle}',
                zorder=(-np.inf if species == 'metric' else None), label=label)
    ax.set_xlim(k_pt[0]/h, k_pt[-1]/h)
    ax.set_ylim(bottom=min(
        np.min(-δρ['g'      ](z)),
        np.min(-δρ['ncdm[0]'](z)),
    ))
    ax.set_xlabel(rf'$k\, \left[h\, \mathrm{{{unit_length}}}^{{-1}}\right]$')
    ax.set_ylabel(
        rf'$-\delta\rho_{{\alpha}}\, '
        rf'\left[{unit_mass}\, \mathrm{{{unit_length}}}^{{-3}}\right]$'
    )
    ax.set_title(rf'$\delta\rho$ perturbations ({gauge} gauge)')
    ax.text(0.8, 0.73, rf'$z = {z}$', transform=ax.transAxes)
axes[1, 0].legend(fontsize=8)
axes[1, 1].yaxis.tick_right()
axes[1, 1].yaxis.set_label_position('right')

# Save figure
fig.tight_layout()
fig.savefig('plot.png', dpi=150)

To run the above script, save it to a file — e.g. plot.py — within the same directory as class_processed.hdf5, then execute it with

python3 plot.py

where python3 is your Python executable.

Tip

The plotting script requires the Python packages NumPy, SciPy, Matplotlib and H5Py, all of which are installed with the Python distribution that comes with CONCEPT. If you have CONCEPT installed (i.e. are not running through Docker), you may then use

source /path/to/concept-installation/concept
$python plot.py

The result of running the script is an image file plot.png with \(2\times 2\) panels, showing

Top left: The evolution history of the background densities \(\bar{\rho}_{\alpha}(a)\) for the different species \(\alpha\), relative to the critical density \(\rho_{\text{crit}}(a)\).
Top right: The matter power spectrum \(P_{\text{b}+\text{cdm}}(k)\) at different redshifts \(z\). This is constructed as

\[P_{\text{b}+\text{cdm}}(z, k) = \zeta^2(k) \delta^2_{\text{b}+\text{cdm}}(z, k)\,,\]

with \(\delta_{\text{b}+\text{cdm}}(z, k)\) obtained from class_processed.hdf5 and \(\zeta(k)\) being the primordial curvature perturbation, parametrised by the primordial parameters \(A_{\text{s}}\), \(n_{\text{s}}\), \(\alpha_{\text{s}}\) and \(k_{\text{pivot}}\) as describe here. The script assigns some values to these parameters. Importantly, these values are independent of the CLASS computation of the \(\delta\) perturbations.
Lower panels: The \(\delta\rho_{\alpha}(k)\) perturbations for the different species \(\alpha\) at two redshifts (early and late). The fictitious \(\delta\rho_{\text{metric}}\) perturbation — supplying general relativistic corrections — generally oscillates wildly, having both positive and negative values. Dotted lines depict negative values that have been flipped.

Note

In CLASS, the sign convention is such that e.g. \(\delta_{\text{cdm}}\) is in fact purely negative. What is plotted is thus \(-\delta_{\alpha}\).

Running through Docker 

As described under Supported platforms, it is possible to run CONCEPT through Docker. While not preferable for running large simulations, running CONCEPT this way for casual use of the class utility is very much encouraged.

Given a parameter file named class-param (e.g. the one above), you may run the CONCEPT class utility through Docker by executing something like

docker run \
    --rm \
    -v "${PWD}:/mnt" \
    jmddk/concept \
        concept \
            -u class \
                --kmin 1e-3/Mpc --kmax 0.1/Mpc --modes 128 \
                --times 512 \
                "b+cdm:1, g" \
            -p /mnt/class-param \
            -n 4

from within the directory of class-param. Given that output_dirs = param.dir is specified within class-param, the resulting HDF5 file will be placed in your current directory.

Note

If running Windows, the above command is valid within PowerShell, not CMD

You can even get rid of the parameter file — running the CONCEPT class utility using a stand-alone one-liner — by supplying the necessary parameters via the -c command-line option.

Tips and tricks 

Below we list a few tips and tricks that are not discussed above.

Plotting perturbations on the fly

It can be difficult to judge whether the various CLASS precision parameters are cranked up high enough to ensure reasonable convergence of the different perturbations. To help with this assessment, it can be useful to see the perturbations plotted, with jagged lines indicating lack of convergence.

By setting the class_plot_perturbations parameter to True, CONCEPT will plot all perturbations entering into the computation of the data in the final class_processed.hdf5 file. E.g.

./concept \
    -u class \
        --kmin <kmin> --kmax <kmax> \
        <perturbations> \    -c "class_plot_perturbations = True"

This will create a class_perturbations and a class_perturbations_processed directory, containing plots of the various perturbations. For details on their content, see the class_plot_perturbations parameter.

The plots within class_plot_perturbations show detrended perturbations, allowing for a detailed view. In these plots, wildly oscillating behaviour is often seen, which may or may not be physical. Non-smooth, jagged behaviour is typically a sign of a lack of convergence, though note that this may not be an issue, since often, most of the information is contained just within the trend-line. To check whether a convergence problem clearly manifest in the full (non-detrended) perturbations, check the plots inside the class_perturbations_processed directory.

Replacing bad CLASS perturbations

Obtaining converged CLASS perturbations for all \(k\) of interest may require extremely high precision CLASS settings, which in turn increases the time and memory requirements for the CLASS computation. Often, converged perturbations in all species are only needed at large scales, as the small scales are completely dominated by matter. Allowing some non-convergence in the other species at high \(k\) is thus often acceptable, with the benefit that lower CLASS precision settings can be used.

However, pathological behaviour in any perturbation may cause trouble (such can be detected using the class_plot_perturbations parameter as described above). If a perturbation is thought to be ill-behaved beyond some \(k\), CONCEPT can be told to replace perturbation values at larger \(k\) with new values obtained via extrapolation from the perturbation values at lower \(k\), using the class_k_max parameter. To e.g. only use the actual CLASS data for \(\partial_{\tau} H_{\text{T}}\) (the conformal time derivative of the trace-free component of the spatial N-body gauge metric) up until \(k = 10\, h\, \text{Mpc}^{-1}\), with data beyond this \(k\) constructed through extrapolation of lower \(k\) data:

./concept \
    -u class \
        --kmin <kmin> --kmax <kmax> \
        <perturbations> \    -c "class_k_max = {'H_T_prime': 10*h/Mpc}"

Ignoring or cleaning the CLASS cache

CONCEPT automatically cashes CLASS computations to disk. When running the CONCEPT class utility several times over with the same input, the first CLASS computation will then be reused, saving computation time.

This behaviour is not always what you want, e.g. if changes to the CLASS source code has been introduced after a disk caching. To request a rerun of CLASS regardless of whether a matching cache is found on disk, set the class_reuse parameter to False, e.g.

./concept \
    -u class \
        --kmin <kmin> --kmax <kmax> \
        <perturbations> \    -c "class_reuse = False"

Even with class_reuse = False, new CLASS computations will be cached to disk. If the disk cache already exists, this will not be overwritten.

To clear the reusable CLASS cache completely, you may run

(source concept && make clean-class-reusable)

class utility

Basic usage

Command-line options

Help: -h, --help

Minimum Fourier mode: --kmin

Maximum Fourier mode: --kmax

Fourier modes: --modes

Times: --times

Gauge: --gauge

File contents

Data layout

Extra output

Units

Complete example (massive neutrinos)

Running through Docker

Tips and tricks

Plotting perturbations on the fly

Replacing bad CLASS perturbations

Ignoring or cleaning the CLASS cache

Basic usage 

Command-line options 

Help: `-h`, `--help`

Minimum Fourier mode: `--kmin`

Maximum Fourier mode: `--kmax`

Fourier modes: `--modes`

Times: `--times`

Gauge: `--gauge`

File contents 

Running through Docker 

Tips and tricks 